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ABSTRACT 


Recent  expansion  in  availability  of  stand-alone  atmospheric  observing  sensors 
introduces  the  question  of  where  placement  maximizes  gain  in  forecast  accuracy. 
This  study  examined  how  sensitivity  analysis  and  observation  targeting  can  be 
used  to  optimize  sensor  placement.  The  primary  objective  of  this  project  was  to 
determine  whether  a  mesoscale  Ensemble  Sensitivity  Analysis  (ESA)  can  be 
used  to  identify  the  sensitivity  profile  of  fog  formation  in  a  complex  terrain 
environment.  Building  on  work  by  Chilcoat  (2012),  this  study  utilized  several 
alternate  methodologies  to  conduct  ESA,  including  a  more  realistic  observing 
network  and  implementation  of  a  Gaussian  filter.  The  second  objective  was  to 
determine  whether  the  calculated  sensitivities  could  be  used  to  reduce  forecast 
uncertainty  for  a  forecast  metric  related  to  dense  fog  formation.  This  was  done 
by  introducing  a  real-world  “truth”  observation  at  the  location  of  greatest 
sensitivity.  The  results  of  this  study  indicate  that  ESA  provides  a  cogent 
mesoscale  sensitivity  profile  which  can  be  used  to  accurately  predict  forecast 
changes  in  fog  using  initial  condition  potential  temperature  values.  This  type  of 
information  may  prove  valuable  in  planning  the  layout  of  future  observational 
networks,  as  well  as  introduce  the  potential  for  performing  data-thinning  during 
data  assimilation  and  real-time  updates  to  forecast  metrics. 
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I.  INTRODUCTION 


A.  MOTIVATION 

The  ability  of  a  weather  forecaster  to  accurately  predict  beyond  the  near 
term  is  predicated  on  the  ability  to  draw  upon  accurate  and  timely  numerical 
weather  prediction  (NWP)  information  to  make  an  informed  meteorological 
outlook.  While  notable  improvements  in  short-term  (less  than  48  hours), 
medium-term  (usually  48-96  hours),  and  long-term  (greater  than  96  hours) 
forecasts  were  made  from  the  late  1970s  to  the  early  1990s,  recent  improvement 
in  short-term  forecasting  has  been  more  limited  (Kalnay  2003).  It  is  believed  that 
advancement  in  NWP  will  likely  come  from  a  variety  of  improvements,  including 
the  creation  of  more  sophisticated  data  assimilation,  improved  ensemble 
forecasting,  employment  of  fully  coupled  atmospheric-hydrological  systems,  and 
more  detailed  atmosphere-ocean-land  coupled  models  (Kalnay  2003).  One 
additional  field  of  improvement  that  is  the  focus  of  this  study  is  the  development 
of  adaptive  observing  strategies,  by  which  “additional  observations  are  placed 
where  ensembles  indicate  there  is  rapid  error  growth  (low  predictability)”  (Kalnay 
2003). 

B.  PROBLEM  STATEMENT 

United  States  Air  Force  commanders,  much  like  their  counterparts  in  the 
private  sector,  are  commonly  constrained  by  fiscal  limitations  when  tasked  to 
execute  operational  objectives.  For  Air  Force  weather  commanders,  these 
objectives  include  the  requirement  to  observe  and  characterize  current  and  future 
weather  conditions  for  military  operations.  The  accuracy  of  these  operational 
products  is  contingent  upon  the  deployment  of  accurate,  reliable  weather  sensors 
as  well  as  access  to  accurate  and  reliable  NWP  products.  These  elements  are 
not  unrelated,  as  an  accurate  meteorological  observation  can  provide  vital 
knowledge  to  a  weather  forecaster  in  the  field  while  also  informing  a  NWP 
system  as  to  the  true  atmospheric  conditions  from  which  to  create  a  more 
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accurate  analysis.  However,  in  a  fiscally  constrained  environment  the  ability  to 
field  and  maintain  the  necessary  number  of  weather  sensors  can  be  impacted, 
thus  limiting  the  flow  of  information  to  both  human  and  computational  forecast 
systems. 

Recent  years  have  seen  the  development  of  a  variety  of  less  expensive 
self-sustaining  weather  sensors.  In  the  U.S.  Air  Force,  a  commonly  used  system 
is  the  “Weather  Pod”,  which  measures  standard  surface  data  (temperature,  dew 
point,  surface  pressure,  wind  speed  and  direction)  and  comes  with  optional 
visibility  and  cloud  sensors.  These  systems  can  be  deployed  to  the  field  at  a 
fraction  of  the  price  of  traditional  tactical  observing  systems,  which  can  exceed 
$100,000. 


Figure  1 .  A  special  operations  weather  team  operator  out  of  Hurlburt  Field 
enables  a  remote  cloud  pod  during  a  training  operation  for  Red  Flag- 
Alaska  10-2  April  15,  2010,  Eielson  Air  Force  Base,  Alaska.  The 
cloud  pod  is  one  of  several  sensors  that  make  up  the  weather  pod. 

(U.S.  Air  Force  photo/Airman  1st  Class  Willard  E.  Grande  II). 

(From  http://www.afweather.af.mil) 

These  self-sustaining  systems  provide  commanders  additional  capability 
to  observe  weather  conditions  at  locations  previously  considered  unobservable, 
either  due  to  enemy  activity  or  lack  of  field  support  (i.e.,  power  or 
communications)  for  observation  systems.  The  positioning  of  these  sensors 
immediately  becomes  a  concern  for  Air  Force  weather  commanders  who  must 
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decide  where  to  place  them  to  maximize  their  usefulness  in  supporting  military 
operations.  In  today’s  battlespace  environment,  these  sensors  likely  accompany 
traditional  atmospheric  sensors,  such  as  the  TMQ-53  Tactical  Meteorological 
Observing  System  (TMOS)  that  may  or  may  not  already  be  deployed  through  the 
rest  of  the  area  of  responsibility  (AOR).  Other  platforms  may  also  be  providing 
atmospheric  data,  including  tactical  weather  radar  (TWR),  remotely  piloted 
aircraft  (RPA)  such  as  the  Predator  vehicle,  and  the  Aerostat,  a  tethered 
intelligence  and  surveillance  airborne  platform  now  common  in  skies  over 
forward  operating  locations.  Placement  of  the  traditional  TMOS  systems  has  in 
the  past  been  determined  by  the  practical  concerns  noted  above.  However,  if 
sensors,  traditional  or  a  more  modern  variety  could  be  placed  in  the  field  that 
would  provide  both  forecasters  and  forecast  systems  critical  data  that  would 
improve  their  accuracy,  this  could  greatly  enhance  meteorological  support  of 
military  operations. 

This  study  examines  the  potential  for  ensemble  sensitivity  analysis  (ESA) 
to  provide  the  capability  to  identify  those  geographic  locations  where  a  more 
accurate  observation  could  most  effectively  improve  a  forecast.  The  Great  Salt 
Lake  (GSL)  Basin  provides  an  ideal  environment  to  test  ESA.  With  its  highly 
varied  terrain,  partially  limited  upstream  observing  network — due  to  the  sparse 
population  density  of  the  Intermountain  West — and  significant  intra-annual 
seasonal  variation,  the  GSL  basin  mimics  a  number  of  operating  environments 
that  the  U.S.  military  operates  in,  including  Afghanistan.  By  proving  the  efficacy 
of  ESA  to  identify  locations  that  provide  the  greatest  impact  on  model  accuracy, 
this  system  may  serve  as  a  model  from  which  future  sensor  networks  may  be 
based  upon  to  ensure  observations  are  taken  at  locations  that  provide  the  most 
useful  atmospheric  information. 

While  ESA  has  shown  itself  to  be  effective  at  the  synoptic  scale,  at  long 
time  scales,  and  in  low  model  resolution,  its  usefulness  in  areas  of  complex 
terrain,  in  shorter  time  scales,  and  at  high  model  resolution  is  less  proven.  One 
of  the  most  important  parameters  to  consider  when  supporting  military  aviation  is 
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visibility,  which  can  significantly  degrade  flying  operations  when  restricted.  The 
most  common  way  for  this  to  occur  is  in  the  presence  of  dense  fog,  which  can 
reduce  visibility  below  mandatory  minimum  operating  thresholds,  typically  1/2 
statute  mile  (SM)  or  less,  depending  on  the  airfield,  the  aircraft  type  and  the 
experience  of  the  pilot  in  command.  In  mountainous  terrain  in  wintertime, 
visibility  is  routinely  reduced  below  operating  threshold  due  to  dense  fog.  The 
formation  and  advection  of  dense  fog  in  the  GSL  Basin  has  been  the  focus  of 
several  studies  (Hill  1987;  Slemmer  2004)  that  identified  the  significant  role  that 
stable  environment  conditions  and  moisture  from  the  GSL  play  in  its  formation. 

This  study  seeks  to  identify  those  locations  that  most  thoroughly  capture 
the  regional  sensitivities  to  stability  and  moisture  exhibited  during  a  dense  fog 
event  in  January  2009.  This  sensitivity  analysis  may  yield  a  more  thorough 
understanding  of  the  predictability  of  fog  forecasting  in  complex  terrain  through 
which  atmospheric  sensors  may  be  placed  more  intelligently,  potentially  allowing 
more  accurate  fog  forecasts  in  similar  environments  throughout  the  world. 

C.  RESEARCH  OBJECTIVES 

The  primary  objective  of  this  project  is  to  further  determine  whether  a 
mesoscale  ESA  can  be  used  to  identify  the  sensitivity  profile  of  fog  formation  in  a 
complex  terrain  environment.  While  this  was  also  the  primary  objective  of 
Chilcoat  (2012),  this  study  utilizes  several  alternate  methodologies  to  achieve  a 
similar  objective,  including  a  more  realistic  observing  network  and  the 
implementation  of  a  Gaussian  filter,  which  is  a  prerequisite  to  properly  utilize  the 
linear  statistics  of  the  ESA  system  of  Torn  and  Hakim  (2008).  The  second 
objective,  also  similar  to  Chilcoat,  is  to  determine  whether  the  calculated 
sensitivities  can  be  used  to  reduce  forecast  uncertainty  for  a  forecast  metric 
related  to  dense  fog  formation.  This  will  be  done  by  introducing  a  real-world 
“truth”  observation  at  the  location(s)  in  the  GSL  basin  with  greatest  sensitivity  to 
determine  if  a  forecast  improvement  can  been  made.  This  type  of  information 
would  prove  highly  valuable  in  planning  the  layout  of  a  future  observational 
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network,  as  well  as  introduce  the  potential  for  performing  data-thinning  during 
data  assimilation  and  real-time  updates  to  forecast  metrics  (Torn  and  Hakim 
2008). 

The  establishment  of  the  viability  of  the  ESA  within  this  mesoscale 
complex  environment  will  be  contingent  on  successful: 

1.  Evaluation  of  the  linearity  of  the  input  variable  (initial  condition)  and 
forecast  variable  through  Gaussianity  testing.  The  assumption  of  linearity  is  a 
prerequisite  to  achieving  accurate  and  reliable  results. 

2.  Introduction  of  a  hypothetical  synthetic  observation  at  the  location  of 
greatest  sensitivity  and  observation  of  a  positive  forecasting  change  through 
increased  forecast  accuracy  (i.e.,  reduced  forecast  error). 

3.  Determination  of  the  ensemble  size  needed  for  ESA  to  provide  accurate 
sensitivity  information  for  this  region. 

Based  upon  the  results  of  this  study,  and  those  presented  by  Chilcoat 
(2012),  a  better  assessment  can  be  made  on  the  usefulness  of  ESA  as  a  tool  for 
improving  prediction  in  complex  terrain  environments. 
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II.  BACKGROUND 


The  world  meteorological  forecasting  community  has  long  relied  on 
numerical  weather  prediction  (NWP)  to  provide  the  primary  means  of  delivering 
accurate  and  reliable  near-to-medium  term  weather  outlooks.  Improvements  in 
NWP  in  recent  decades  can  be  attributed  to  the  increased  power  of 
supercomputers,  the  improved  representation  of  small-scale  physical  processes, 
the  use  of  more  accurate  methods  of  data  assimilation  (DA),  and  increased 
availability  of  data  (Kalnay  2003).  The  improvement  in  data  assimilation  methods 
has  enabled  a  more  accurate  representation  of  the  real  atmosphere,  thus 
providing  more  realistic  initial  conditions  for  model  integration,  and  more  accurate 
forecasts. 

A.  DATA  ASSIMILATION  AND  THE  DEVELOPMENT  OF  THE  ENSEMBLE 

KALMAN  FILTER 

The  importance  of  DA  for  the  accuracy  of  NWP  has  been  established 
since  at  least  the  1950s.  NWP  is  an  initial-value  problem  that  considers  a 
representation  of  the  real  atmosphere  and  integrates  it  forward  to  simulate  the 
real  atmosphere  (Kalnay  2003).  Because  it  is  necessary  to  characterize  the 
atmosphere  at  significantly  more  data  points  than  can  be  realistically  observed,  a 
methodology  was  developed  to  “spread”  the  known  observational  information  to 
the  unobserved  regions.  Later,  this  process  came  to  further  include  the 
correction  of  the  model  estimate  of  the  atmosphere  when  it  differed  from  the 
observation.  This  process,  also  known  as  performing  the  analysis,  became  the 
foundation  for  improving  DA  in  the  1950s  and  1960s.  Early  methods  to  perform 
this  function  included  the  use  of  climatology,  the  successive  correction  method 
(SCM)  and  optimal  interpolation  (Ol)  (Kalnay  2003). 

A  significant  step  forward  came  with  the  introduction  of  the  three- 
dimensional  variational  data  assimilation  (3DVAR)  system  in  the  early  1990s 
(Kalnay  2003).  This  system,  based  largely  on  the  work  of  Sasaki  from  the  late 
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1950s  through  early  1970s  (Lewis  and  Lakshmivarahan  2008),  used 
understanding  of  established  dynamical  relationships  to  identify  covariance 
relationships  which  could  in  turn  be  used  to  perform  the  model  analysis.  These 
dynamical  relationships  were  determined  through  sampling  methods  over  long 
periods  of  time  and  utilizing  multiple  case  studies  (Hakim  and  Torn  2008).  The 
limitation  of  this  system  lies  in  its  deterministic  method  of  performing  state 
estimation,  which  fails  to  account  for  known  errors  in  initial  conditions  that  have 
identifiable  limitations  on  predictability  (Lorenz  1963). 

An  important  advancement  was  made  when  Evensen  (1994)  successfully 
implemented  a  Kalman  filter  using  ensemble  techniques  to  determine  the 
forecast-error  statistics.  This  technique  is  notably  different  from  3DVAR  and  is 
known  as  an  ensemble  Kalman  filter  (EnKF).  The  EnKF  provides  a  systematic, 
temporally  evolving  means  of  determining  the  error  statistics,  which  are 
essentially  the  forecast-error  statistical  covariances.  This  data  may  then  be  used 
to  adaptively  alter  the  DA  process  to  provide  greater  emphasis  (more  weighting) 
to  those  model  areas  with  less  statistical  error,  and  less  emphasis  (less 
weighting)  to  those  model  areas  with  more  statistical  error.  By  providing 
probabilistic,  flow-dependent  estimates  of  error  statistics,  the  EnKF  provides  the 
basis  for  the  technique  used  in  this  study,  ensemble  sensitivity  analysis  (ESA). 

B.  SENSITIVITY  ANALYSIS 

While  success  was  being  made  in  the  implementation  of  these  ensemble 
methods  in  the  late  1990s,  and  early  2000s,  the  area  of  sensitivity  analysis  was 
attracting  considerable  research.  By  definition,  this  involves  the  study  of  how 
uncertainty  in  model  output  (numerical  or  otherwise)  can  be  related  to  differing 
sources  of  uncertainty  in  the  model  input  (Saltelli  et  al.  2004).  Within  a 
meteorological  context,  this  began  by  seeking  to  understand  how  chosen  aspects 
of  a  synoptic  feature  may  change  due  to  alterations  in  physical  or  dynamical 
processes,  or  initial  atmospheric  conditions  (Errico  1997).  Analysis  of  this  nature 
offered  the  opportunity  to  provide  insight  into  relationships  between  model  input 
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and  output,  possibly  identifying  model  forecast  error  that  may  be  caused  by 
boundary  or  initial  conditions.  A  secondary  use  was  to  determine  what  initial 
synoptic  features  most  strongly  related  to  dependent  synoptic  features  are  later 
times  (Langland  et  al.  1995).  Identifying  these  relationships  may  help  to  provide 
a  better  understanding  of  the  dynamics  of  forecast  error,  as  well  as  identify 
possible  locations  to  place  additional  observations  to  reduce  forecast  error  (Torn 
and  Hakim  2008). 

Initial  advancement  in  this  field  focused  on  measuring  changes  in  forecast 
parameters  by  manipulating  single  input  variables.  This  method,  though  simple 
to  implement,  failed  to  accurately  identify  the  impact  of  all  possible  input 
perturbations  on  predicted  forecast  metrics,  and  its  utility  was  determined  to  be 
limited  (Errico  1997).  A  more  advanced  method  was  then  developed  to  perform 
sensitivity  analysis  which  relied  upon  the  adjoints  of  tangent-linear  models.  The 
adjoint  of  a  NWP  model  can  be  thought  of  as  the  transpose  of  the  tangent  linear 
of  the  NWP  model  itself.  However,  the  limitations  and  difficulty  of  developing  and 
maintaining  tangent  linear  models  is  well  documented  (Hakim  and  Torn  2008; 
Torn  and  Hakim  2008;  Ancell  and  Hakim  2007).  Recent  research  has  shown 
promise  in  computing  the  sensitivity  analysis  in  a  more  straightforward  manner. 

This  advancement  involves  the  use  of  ensembles,  rather  than  adjoints,  to 
compute  sensitivity  analysis.  This  technique,  ensemble  sensitivity  analysis 
(ESA),  uses  ensembles  and  probabilistic  rather  than  deterministic  analyses  to 
estimate  relationships  between  forecast  metrics  and  initial  conditions.  This  area 
of  study  has  its  roots  in  targeting,  or  adaptive  observational  experiments,  which 
seek  to  understand  where  analysis  and  forecast  improvement  could  be  achieved 
through  placement  of  additional  observations  (Hamill  and  Snyder  2002).  In  the 
late  1990s  and  early  2000s  a  multitude  of  studies  were  conducted  and  ideas 
proposed  that  sought  the  most  efficient  and  valid  means  of  solving  this  problem. 
These  methods  included  the  singular  vector  technique,  ensemble  spread 
technique,  ensemble  transform  technique,  and  the  ensemble  transform  Kalman 
filter,  among  several  others  (Hamill  and  Snyder  2002).  Building  on  the  concepts 
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put  forth  by  Berliner  et  al.  (1999)  regarding  forecast  and  analysis  error  statistics, 
and  by  Burgers  et  al.  (1998)  regarding  background  error  covariance  techniques, 
Hamill  and  Snyder  proposed  their  own  theoretical  method  of  targeting 
observations  using  an  EnKF  to  minimize  analysis  errors  (Hamill  and  Snyder 
2002). 

C.  ENSEMBLE  SENSITIVITY  ANALYSIS 

In  2005,  a  methodology  based  upon  the  EnKF  and  the  sensitivity  analysis 
work  described  above  was  codified  as  ESA,  and  was  demonstrated  to  be  a 
viable  means  of  calculating  state  estimation  (e.g.,  expected  value),  state 
uncertainty  (e.g.,  variance)  and  the  relationship  between  all  variables  at  all 
locations  in  time  (Hakim  and  Torn  2008).  In  their  proposal,  Hakim  and  Torn 
devised  a  method  to  compute  the  sensitivity  of  a  forecast  metric  to  changes  in 
model  initial  conditions  using  Gaussian  statistics  and  linear  regression  of  model 
statistics  determined  by  the  ensemble  members.  In  ESA,  the  uncertainty 
identified  by  Lorenz  (1963)  that  accompanies  the  propagation  of  imperfect 
observations  under  nonlinear  chaotic  dynamics  is  accounted  for  by  the 
production  of  instantaneous  and  independent  state  estimations  that  are  provided 
by  the  ensemble  analysis  (Gombos  and  Hansen  2008).  The  utility  of  this  method 
is  that  it  provides  a  means  of  determining  the  possible  implications  of  forecast 
error,  and  where  this  error  may  be  reduced  through  active  management  of 
observation  locations  (Torn  and  Hakim  2008).  In  other  words,  this  technique 
provides  an  optimal  means  of  performing  observation  targeting.  This  technique 
has  shown  efficacy  as  an  alternative  to  adjoint  methods  (Figure  2,  Ancell  and 
Hakim  2007),  and  also  a  successful  means  of  determining  forecast  error  in 
synoptic  storms  (Hakim  and  Torn  2008),  storms  undergoing  extratropical 
transition  (Torn  and  Hakim  2009),  and  in  developing  hurricanes  (Torn  2010). 
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Figure  2.  Plot  showing  the  effective  equivalence  of  adjoint  and  ensemble 

sensitivity  methodologies.  Plotted  is  the  change  in  response  function 
associated  with  both  the  statistical  projection  onto  the  adjoint  sensitivity  field 
(denoted  by  the  symbol  o)  and  the  ensemble  sensitivity  (denoted  by  the  symbol 
x)  versus  the  nonlinear  change  in  response  function  for  a  statistically  spread  unit 
temperature  perturbation  made  at  20  grid  points  throughout  the  model  domain. 

(From  Ancell  and  Flakim  2007) 

Though  successfully  demonstrated  as  a  viable  alternative  to  adjoint-based 
methods,  these  studies  applied  ESA  only  on  the  synoptic  scale  and  with 
relatively  coarse  models.  Recent  studies  have  begun  to  examine  the  efficacy  of 
ESA  for  use  in  finer  resolution  models  and  at  smaller  horizontal  scales 
(mesoscale  and  below).  Chilcoat  (2012)  applied  ESA  techniques  using  an  EnKF 
and  high-resolution  forecast  models  (4-km  grid  spacing)  in  complex  terrain  (Great 
Salt  Lake  Basin).  Chilcoat  demonstrated  ESA  as  a  useful  methodology  in 
accurately  identifying  geographic  locations  where  sensitivity  between  an  initial 
state  variable  (u-wind)  and  a  forecast  parameter  (fog,  represented  as  water 
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vapor  mixing  ration,  qv)  was  greatest.  The  next  step,  and  the  purpose  of  this 
research,  is  to  validate  the  sensitivities  identified  by  Chilcoat  and  then  to 
determine  if  a  more  accurate  forecast  can  be  created  by  adding  an  observation 
at  these  locations. 


12 


III.  DATA  AND  METHODOLOGY 


The  methodology  used  in  this  research  is  based  upon  the  ensemble 
sensitivity  technique  defined  in  the  Background  section  of  this  thesis.  To  fully 
understand  how  ESA  works  and  how  it  was  used  for  this  research,  it  is  necessary 
to  first  identify  the  procedure.  Gombos  and  Hansen  describe  ESA  as  a  method 
of  “statistically  inferring  synoptic  dynamical  relationships  from  the  covariability  of 
fields’  ensemble  analysis”  (Gombos  and  Hansen  2008).  Stated  another  way,  this 
method  seeks  to  identify  statistical  relationships  between  observations  (initial 
conditions)  and  forecast  fields  and  then  to  draw  conclusions  from  the  results  that 
can  be  equated  to  the  dynamic  relationships  that  exist  between  the  fields.  This 
technique  is  part  of  an  emerging  field  of  study  that  uses  underlying  error  statistics 
such  as  the  variance  and  covariance  of  error  in  observations,  model  analysis, 
and  model  background  fields,  to  identify  relationships  within  complex  nonlinear 
systems.  These  statistical  relationships  are  based  upon  linear  regression  and 
therefore  are  subject  to  the  limitations  of  linear  statistics  which  assume  normally 
distributed  data  fields,  with  variation  in  the  distributions  attributable  to  normally 
distributed  random  error. 

A.  LINEARITY  TESTING 

To  test  these  assumptions  and  thus  the  validity  of  the  linear  regressions 
used  in  this  technique,  the  Gaussianity  of  the  distributions  must  be  examined. 
Wilks  defined  linear  regression  as  an  attempt  to  “summarize  the  relationship 
between  two  variables,  shown  graphically  in  their  scatterplot,  by  a  single  straight 
line”  (Wilks  2006).  Simple  linear  regression  defines  the  linear  relationship  that 
exists  between  two  variables,  typically  x  (independent  or  predictor  variable)  and  y 
(dependent  or  predictand  variable).  This  linear  relationship  is  defined  by  the 
equation  y  =  mx  +  b  where  m  is  the  slope  of  the  line  and  b  is  the  y-intercept.  The 
line  that  is  created  during  this  procedure,  the  best  fit  line,  will  produce  the  least 
error  in  expected  value  of  y,  when  given  values  of  x.  This  error  minimization  is 
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achieved  by  finding  the  lowest  sum  of  squared  errors  between  the  values  of  y, 
after  inputting  the  values  for  x.  This  is  also  known  as  ordinary  least  squared 
regression  (OLS). 

The  assumption  made  with  linear  regression  is  that  the  given  x  and  y 
variables  are  normally  distributed,  also  known  as  a  Gaussian  distribution.  To 
determine  whether  or  not  these  variables  are  Gaussian,  we  applied  a  Lilliefors 
test  to  the  variable  sets.  The  Lilliefors  test,  a  variation  of  the  Kolmogrov-Smirnoff 
test,  is  a  two-sided  goodness-of-fit  test  that  measures  the  distribution  of  each 
variable  set  (empirical)  to  determine  how  well  it  compares  against  the  known 
normal  distribution  (theoretical).  To  perform  this  test,  we  propose  the  null 
hypothesis  that  the  empirical  data  comes  from  a  normally  distributed  population. 
The  alternative  is  that  the  data  does  not  come  from  a  normal  distribution.  The 
test  statistic  used  is  defined  as: 

D„  =  max  |Fn(x)-F(x)| 

where  Fii(x)  is  the  empirical  cumulative  distribution  function  (CDF)  estimated 
from  the  sample  and  F(x)  is  the  normal  CDF.  The  mean  and  standard 

deviation  of  both  CDFs  must  be  equal.  The  test  statistic  Dn  determines  the 

largest  total  difference  between  the  empirical  and  theoretical  CDF.  The 
computed  test  statistic  is  then  related  to  a  critical  value  using  a  table.  Matrix 
Laboratory  (MATLAB)  uses  a  table  of  critical  values  that  have  been  computed 
using  a  Monte  Carlo  simulation  for  sample  sizes  less  than  1000,  and  significance 
levels  (a)  between  0.001  and  0.50.  If  after  performing  the  test  the  total  difference 
between  the  empirical  and  theoretical  CDF’s  is  too  great,  a  test  statistic  is 
returned  that  exceeds  the  significance  level  and  the  null  hypothesis  will  be 
rejected.  When  this  occurs,  the  data  is  rejected  for  its  failure  to  be  reliably  drawn 
from  a  normally  distributed  population.  The  significance  level  chosen  for  this 
study  for  this  test  is  0.05. 
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The  Lilliefors  test  in  this  study  was  completed  upon  both  the  observed 
variable  (x)  and  the  forecast  metric  (J),  as  the  assumption  of  Gaussianity  applies 
equally  to  both  variables.  Non-Gaussian  observed  and  forecast  metrics  were 
rejected  if  they  failed  to  pass  the  the  Lilliefors  test.  By  selecting  only  data  sets 
that  have  Gaussian  distributions,  we  hoped  to  ensure  the  assumptions  of  linear 
statistics  are  upheld. 

B.  ESA  CALCULATION 

The  linear  regression  defined  in  3.A.  was  identified  by  Torn  and  Hakim  in 
2008  as  ESA,  a  calculation  of  the  sensitivity  of  the  forecast  metric  to  the  initial 
conditions  (Torn  and  Hakim  2008).  The  authors  showed  that  sensitivity  of  a 
forecast  metric  J  to  an  initial  analysis  state  variable  xt  is  defined,  in  scalar  form, 
by 

dJ  _  cov(y.x) 
dxi  var(x)(. 

Covariance,  which  quantifies  the  strength  of  the  linear  relationship 
between  the  two  arguments,  is  denoted  by  cov  while  var,  which  quantifies  how 
far  the  population  is  spread  from  the  mean,  is  the  variance.  In  scalar  form,  x( 
represents  the  analysis  state  variable  for  a  single  grid  point.  Any  physical 
variable  can  be  represented  by  xi ,  samples  of  which  are  provided  by  an 

ensemble  of  size  M.  Provided  that  an  EnKF  is  used  for  assimilating 
observations,  this  ensemble  system  provides  estimates  of  the  ensemble  mean 
values  of  state  variable  and  forecast  metric,  from  which  sensitivities  across  the 
entire  grid  may  be  computed. 

To  further  identify  the  linear  dependence  between  the  two  variables  J  and 
x,  the  Pearson  product-moment  coefficient  of  linear  correlation,  or  correlation 

coefficient  (Wilks  2006),  is  introduced  to  the  ESA  equation.  This  value  is  defined 
by  the  equation 
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cov(/,x) 
r  — - — 

where  r,  a  unit-less  measure,  is  the  correlation  coefficient  between  J  and  xn 
while  a,  and  <jx  represent  the  standard  deviation  of  the  sample  of  the 
respective  variables.  This  equation  can  be  rewritten  as 

cov(/,  x  )  =  rJ  x  <7j  (<ta  . 

If  we  then  substitute  this  equation  into  our  original  ESA  equation,  the  updated 
form,  with  var(x),  rewritten  in  symbolic  form,  appears  as 


CJ  _r,aJ(cTx)i 


This  after  cancellation  becomes 

dj  _  r^Qj 

fa,  K), 

This  new  form  of  Hakim  and  Torn  (2008)  can  be  used  to  apply  a  confidence 
interval  test  to  account  for  sampling  error. 

This  study  seeks  to  define  the  sensitivity  dJ/dx,  between  the  observed 
variable  x,-  represented  as  a  column  state  vector  x  containing  near  surface 
temperature  values  and  the  forecast  metric  J  as  the  water  vapor  mixing  ratio,  qv, 
at  our  point  of  interest,  the  Salt  Lake  City  International  Airport. 

C.  STATISTICAL  SAMPLING  ERROR 

Because  the  ensemble  system  draws  upon  a  finite  sample,  sampling  error 
is  naturally  introduced  into  the  ESA.  As  identified  by  Hakim  and  Torn  (2008),  it  is 
necessary  to  ensure  that  the  determined  sensitivities  are  statistically  significant  in 
order  to  limit  sampling  error. 

To  account  for  this  sampling  error,  we  will  perform  a  confidence  test  at  the 
95%  confidence  interval,  for  which  the  null  hypothesis  is  that  changes  to  input 
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variable  x;  do  not  change  the  forecast  metric  J.  We  will  reject  this  hypothesis  if 
the  absolute  value  of  the  correlation  coefficient  identified  in  3.B.  exceeds  the  95% 
confidence  bounds  as  estimated  by  the  ensemble  data  (Torn  and  Hakim  2008). 

We  seek  to  define  the  probability,  or  p-value,  which  reflects  the  measure  of 
evidence  against  our  formulated  null  hypothesis.  A  small  p-value,  p[J,  x],  will 
correspond  to  strong  evidence  that  we  can  accept  the  correlation  coefficient, 
r[J,  Xi]  ,  which  accompanies  each  sensitivity  calculation.  Working  within  MATLAB 
the  p-values  were  determined  by  utilizing  the  correlation  coefficients  to  generate 
a  test  statistic  with  N-2  degrees  of  freedom.  These  associated  values  are  then 
used  to  test  the  null  hypothesis.  With  a  95%  confidence  interval  (Cl),  when  low 
p-values  (less  than  0.05)  are  calculated,  we  can  conclude  that  the  correlation 
coefficients  are  statistically  significant,  which  equates  to  the  rejection  of  the  null 
hypothesis.  This  test  is  performed  for  each  sensitivity  calculation  and  locations 
with  correlation  coefficients  that  fail  the  statistical  significance  test  are  rejected 
and  data  provided  from  these  points  is  masked. 

D.  MODEL  CONFIGURATION 

The  ensemble  system  used  in  this  research  is  the  Advanced  Research 
version  of  the  Weather  Research  and  Forecasting  model  (WRF-ARW,  version 
3.2.1).  The  unique  contribution  of  the  ARW  subset  comes  from  the  addition  of  a 
dynamics  solver  and  specific  physics  schemes,  numerical/dynamics  options, 
initialization  routines,  and  a  data  assimilation  package  (Skamarock  et  al.  2008). 

1.  Nesting  Configuration 

The  ensemble  system  configuration  will  mirror  that  used  by  Chilcoat 
(2012),  with  a  telescoping  fixed  nested  grid  with  36-km  horizontal  grid  spacing  on 
the  outermost  domain,  12-km  horizontal  grid  spacing  on  the  middle  domain,  and 
4-km  grid  spacing  on  the  innermost  domain.  All  domains  in  the  model  contain  60 
vertical  /7-levels,  each  with  upper  atmospheric  boundary  set  at  100  hPa.  The 
telescoping  grid  domains  are  pictured  in  Figure  3.  The  innermost  domain  lies 
over  the  central  Intermountain  West  and  is  centered  slightly  to  the  west  of  the 
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Great  Salt  Lake.  The  dimensions  of  this  inner  domain  include  222  x  183  grid 
points  which  cover  an  area  888  km  x  732  km. 


Figure  3.  Telescoping  nested  configuration  used  in  this  study.  The  outermost 
grid  domain  uses  36-km  horizontal  grid  spacing,  the  middle  domain  uses  12-km 
horizontal  grid  spacing  and  the  innermost  domain  uses  4  km  horizontal  grid 

spacing.  (From  Chilcoat  2012) 

2.  Vertical  Coordinate  System 

The  WRF-ARW  vertical  coordinate  system  has  been  established  with 
terrain-following,  dry  hydrostatic  pressure,  with  vertical  grid  stretching  permitted 
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(Skamarock  et  al.  2008).  The  vertical  coordinates  are  denoted  by  q  and  vary 
from  a  value  of  1  at  the  surface  to  0  at  the  upper  boundary  of  the  model  domain. 
Figure  4  shows  a  visual  representation  of  the  vertical  coordinate  system. 
Because  the  primary  focus  of  this  study  was  dense  fog  formation  and  advection 
at  the  surface  and  near-surface  level,  the  two  lowest  /7-levels  were  of  most 
importance.  Table  1  shows  the  approximate  height  of  the  /7-levels  in  the  lowest 
500  meters  of  the  atmosphere. 


Figure  4.  ARW  q  coordinate  (From  Skamarock  et  al.  2008) 
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n-Level 

AGL  Height  (m) 

1 

25 

2 

60 

3 

107 

4 

166 

5 

240 

6 

334 

7 

446 

8 

513 

Table  1 .  Approximate  above-ground  level  (AGL)  height  for  /7-levels  in  the 

lowest  500  meters  of  the  atmosphere. 

3.  Physics 

WRF-ARW  separates  physics  options  into  five  categories,  each  of  which 
contain  several  options.  The  primary  physics  categories  include  (1) 
microphysics,  (2)  cumulus  parameterization,  (3)  planetary  boundary  layer  (PBL), 
(4)  land-surface  model,  and  (5)  radiation  (Skamarock  et  al.  2008).  For  this  study, 
the  physics  options  selected  are  the  same  as  those  used  operationally  by  the  Air 
Force  Weather  Agency  (AFWA)  in  their  deterministic  WRF  model.  The 
microphysics  scheme  used  for  this  study  is  the  WRF  Single-Moment  5-class 
(WSM5)  scheme.  This  scheme  allows  water  to  be  held  in  five  separate 
hydrometeorological  categories  and  is  efficient  in  intermediate  grids  between  the 
cloud-resolving  grids  and  the  mesoscale  grids  (Skamarock  et  al.  2008).  Also 
used  is  the  Kain-Fritsch  cumulus  parameterization  (Kain  2004),  though  it  was 
enabled  on  the  two  outermost  grids  only,  as  cumulus  parameterizations  can  be 
considered  theoretically  valid  only  on  course  grid  scales,  generally  greater  than 
10  km  (Skamarock  et  al.  2008).  The  planetary  boundary  layer  (PBL)  scheme 
implemented  is  the  Yonsei  University  PBL  scheme  (YSU  PBL)  (Hong  et  al. 
2006),  while  the  land  surface  model  is  the  Noah  (Ek  et  al.  2003).  Finally,  the 
radiation  parameterization  scheme  implemented  was  the  rapid  radiative  transfer 
model  (RRTM;  Mlawer  et  al.  1997)  for  longwave  processes  and  the  Dudhia 

20 


shortwave  scheme  (Dudhia  1989)  for  shortwave  processes.  Table  2  summarizes 
the  physics  parameterizations  selected  for  this  study. 


Microphysics 

Cumulus 

PBL 

LSM 

Longwave 

Shortwave 

WSM5 

Kain-Fritsch 

YSU 

Noah 

RRTM 

Dudhia 

Table  2.  WRF-ARW  physics  options  selected  for  this  study. 


4.  Data  Assimilation  Process 

Data  Assimilation  (DA)  was  performed  with  the  Data  Assimilation 
Research  Testbed  (DART).  DART  is  a  community  software  environment  created 
for  ensemble  DA  research,  developed  and  maintained  at  the  National  Center  for 
Atmospheric  Research  (NCAR)  (Anderson  et  al.  2009).  For  this  project  DART 
was  run  on  Cray  XE6  supercomputers  at  the  U.S.  Army  Engineer  Research  and 
Development  Center  (ERDC)  in  Vicksburg,  Mississippi. 

a.  Truth  Run 

To  recreate  the  conditions  that  supported  this  dense  fog  event,  an 
observing  system  simulation  experiment  (OSSE)  was  conducted.  The  purpose 
of  this  experiment  was  to  provide  a  WRF-based  synthetic  “truth”  representation 
which  resembled  the  true  state  of  the  atmosphere  as  closely  as  possible.  The 
WRF  “truth”  simulation  was  initialized  at  0000  UTC  on  19  January  2009  using 
conditions  provided  by  the  North  American  Regional  Reanalysis  (NARR).  Lateral 
boundary  conditions  were  drawn  from  NARR,  with  perturbations  implemented  per 
Torn  et  al.  (2006).  From  this  point,  the  WRF  soil  conditions  were  reset  every 
three  hours  to  NARR-provided  values  to  ensure  surface-based  boundary  layer 
forcing  mimicked  reality  as  closely  as  possible.  This  simulation  provided  a 
database  to  test  the  results  of  our  experiments  against. 

b.  Ensemble  System 

Conditioning  of  the  WRF  ensemble  system  began  at  0000  UTC  on 

20  January.  DA  proceeded  every  three  hours  through  the  introduction  of 
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available  “truth”  rawindsonde  observations  that  were  recorded  across  the  outer 
grid  domain  during  this  time.  While  Chilcoat  (2012)  introduced  synthetic 
rawindsonde  observations  at  100  random  locations  across  the  domain  to 

simulate  a  sparse  data  environment,  the  approach  in  this  study  provides 

observations  that  spatially  represent  the  rawinsonde  data  network  present  in 
western  North  America.  Besides  the  introduction  of  lateral  boundary  conditions 
and  soil  reset  from  NARR,  no  other  observations  were  added  to  the  system.  For 
two  days,  the  ensemble  was  integrated  forward  until  0000  UTC  on  22  January 
2009,  when  it  was  concluded,  as  shown  by  Torn  et  al.  (2006),  that  the 

perturbations  entering  through  the  lateral  boundaries  were  being  properly 

corrected  within  the  limited  area  domain  of  the  ensemble  system. 

At  initialization  time  the  96-members  of  the  ensemble  system  were 
formed  by  adding  random,  yet  spatially  consistent  perturbations  to  the  truth. 
Perturbations  were  drawn  from  the  static  error  covariance  field  for  the  WRF-Var. 
Every  three  hours  available  “truth”  rawindsonde  observations  were  assimilated 
again,  enabling  the  creation  of  an  entirely  new  96-member  ensemble  system. 
Perturbations  drawn  from  the  WRF-Var  covariance  fields  were  applied  to  the 
lateral  boundary  conditions  to  ensure  that  perturbations  continued  to  enter  the 
limited  area  domain  from  outside. 

Initial  conditions  from  the  DA  were  then  used  to  create  forecasts 
that  extended  past  the  assimilation  cycling  window.  Ensemble  forecasts  were 
created  at  1200  UTC  23  January  and  1800  UTC  23  January.  Both  systems  were 
designed  to  end  at  0000  UTC  24  January,  approximately  when  the  dense  fog 
event  began. 

In  order  to  examine  how  the  size  of  the  ensemble  system  affected 

the  ESA,  reduced-size  ensemble  systems  were  created  alongside  the  full  96- 

member  field.  At  0000  UTC  on  22  January,  sub-samples  with  80,  64,  and  48 

ensemble  members  were  created  and  run  separately,  yet  analogously  to  the  96- 

member  system.  These  sub-systems  were  created  by  selecting  the  first  80,  64, 

and  48  ensemble  members  of  the  full  96-member  ensemble  system.  By  drawing 
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all  ensemble  systems  from  the  same  96-member  background  field,  we  can  be 
certain  that  our  sensitivities  are  drawn  from  the  same  ensemble  data  system. 
This  limits  the  variation  in  the  initial  conditions  of  the  systems  and  provides  a 
greater  likelihood  of  model  systems  with  similar  parcel  trajectories  and  therefore 
similar  sensitivities.  All  four  ensemble  systems  end  at  0000  UTC  on  24  January. 
Figure  5  represents  the  full  DA  cycle. 
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Figure  5.  Data  Assimilation  Cycle  used  for  this  project.  All  forecasts  are  valid 
at  0000  UTC  24  January  2009  to  coincide  with  the  onset  of  the  dense  fog  event 

at  Salt  Lake  City  International  Airport. 

c.  Covariance  Localization  and  Covariance  inflation 

Because  of  sampling  error  inherent  in  the  ensemble  system,  it  is 
necessary  to  implement  a  tool  to  limit  the  horizontal  and  vertical  magnitude  of  the 
covariance  between  x,  and  J.  This  technique  is  covariance  localization  and  helps 
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prevent  the  creation  of  spurious  covariances  at  long  distances  that  come  about 
because  of  the  inherent  model  error.  While  this  process  does  introduce  some 
imbalances  to  the  system,  it  is  a  “practical  necessity  for  EnKF  NWP  applications” 
(Hamill  2013). 

Another  corrective  tool  implemented  was  covariance  inflation.  This 
technique  increases  the  spread  of  the  ensemble  in  order  to  negate  any  possible 
underestimation  of  forecast  error  variance  that  may  occur  because  of  limitations 
in  ensemble  size.  This  correction  helps  to  ensure  that  ensemble  variance  is  not 
artificially  limited  because  of  model  error. 

E.  REGIONAL  OVERVIEW 

1 .  Geography  of  Great  Salt  Lake  Basin 

This  research  focuses  on  the  Great  Salt  Lake  Basin  of  Utah,  with  an 
emphasis  on  conditions  near  the  Great  Salt  Lake  (GSL)  and  at  Salt  Lake  City 
International  Airport  (ICAO:  KSLC).  The  GSL  Basin  is  a  snow-dominated 
hydrologic  system  situated  in  northern  Utah  that  encompasses  approximately 
55,000  sq  km  in  area.  The  GSL  covers  an  average  area  of  4,500  sq  km,  though 
it  can  exhibit  significant  (+/-  1 ,300  sq  km)  interannual  variability  due  to  its  shallow 
depth  and  dependence  on  hydrologic  cycles.  Because  of  its  shallow  depth  and 
high  salinity,  water  temperatures  vary  from  subfreezing  in  the  winter  to  greater 
than  80  °F  in  the  summer.  The  average  elevation  of  the  GSL  is  4,202  feet  above 
sea  level  and  it  is  bounded  by  terrain  in  nearly  all  directions  (Figure  6).  To  the 
east,  and  most  prominently,  is  the  Wasatch  Range  (Figure  6,  Point  1)  which  runs 
north  to  south  along  the  length  of  the  lake  and  ranges  in  elevation  from  4,200  to 
nearly  12,000  feet,  with  most  crests  above  10,000  feet.  Across  the  south,  west, 
and  north  edges  of  the  GSL  are  a  series  of  smaller,  primarily  north  to  south 
oriented  mountains  which  serve  to  further  bound  the  lake  by  creating  a  vertical 
barrier  that  inhibits  interaction  with  the  rest  of  the  intermountain  west.  Within  or 
near  the  southern  boundaries  of  the  GSL,  these  peaks  include  the  Oquirrh 
Mountains  (Figure  6,  Point  2),  which  rises  to  10,620  feet,  Antelope  Island  (Figure 
6,  Point  3)  highest  elevation  6,596  feet,  Stansbury  Island  (Figure  6,  Point  4), 
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highest  elevation  6,647  feet,  and  the  Stansbury  Mountains  (Figure  6,  Point  5) 
further  to  the  south,  which  reaches  an  elevation  of  11,035  feet  at  Deseret  Peak. 
To  the  southwest  of  the  GSL  are  the  Lakeside  Mountains  (Figure  6,  Point  6), 
capped  by  Craner  Peak  at  6,621  feet,  while  to  the  northwest  sits  the  Hogup 
Mountains  (Figure  6,  Point  7),  which  exceeds  7,050  feet  at  Tangent  Peak. 
Finally,  north  of  the  GSL  lies  a  series  of  north  to  south  ranges  which  extend  north 
to  Idaho,  the  most  prominent  of  which  are  the  Promontory  Mountains  (Figure  6, 
Point  8),  rising  over  7,300  feet. 


Figure  6.  Map  of  Great  Salt  Lake  with  primary  terrain  features  annotated. 
(After  Utah  Automated  Geographic  Resource  Center,  http://gis.utah.gov/). 

The  cumulative  effect  of  these  mountains  is  the  creation  of  a  vertical 


barrier  in  excess  of  2,000  feet  around  the  GSL  that  impedes  interaction  with  the 
rest  of  the  Intermountain  West.  The  primarily  north  to  south  orientation  of  these 
orographic  features  tends  to  block  zonal  flow  and  channel  meridional  flow  across 
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the  GSL  basin.  In  the  fall,  spring  and  winter,  large  gaps  between  the  north  to 
south  oriented  features  enable  the  advection  of  cold  air  from  southern  Idaho  and 
warm  advection  from  points  southward,  while  disrupting  west  to  east  oriented 
flows.  In  the  winter  as  the  lake  cools,  the  flow-blocking  provided  by  these 
mountains  helps  enable  the  buildup  of  a  dome  of  cold  air  across  the  GSL  basin, 
which  plays  a  significant  role  in  the  lifecycle  of  mesoscale  meteorological 
phenomena,  including  the  formation  of  deep  stable  layers  (DSL)  as  identified  by 
Wolyn  and  McKee  (1988). 

Salt  Lake  City  International  Airport  (KSLC)  lies  southeast  of  the  GSL 
(Figure  6,  black  star)  approximately  eight  miles  from  boundary  of  the  lake.  The 
airport  is  located  four  miles  northwest  of  downtown  Salt  Lake  City.  The  airfield 
elevation  is  1 ,288  meters  (4,227  feet)  which  situates  KSLC  25  feet  higher  than 
the  GSL,  on  average. 

2.  Fog  at  Salt  Lake  City  International  Airport 

Though  the  general  dynamics  of  fog  formation  in  valley  basins  is  well 
understood,  analysis  of  fog  formation  in  the  GSL  basin  is  not  widely  available. 
Published  papers  on  the  subject  are  limited  and  include  a  study  on  the  effect  of 
the  size  of  the  GSL  on  fog  formation  (Hill  1988)  as  well  as  partial  examination  of 
fog  formation  mechanisms  in  Wolyn  and  McKee’s  examination  of  deep  stable 
layers  (Wolyn  and  McKee  1988).  Unpublished  materials  are  more  numerous  and 
include  a  study  on  dense  fog  initiation  in  the  Salt  Lake  Valley  by  Hogan  (1998) 
and  a  database  of  fog  events  by  Alder,  which  was  included  in  a  broader 
climatological  publication  on  historical  weather  records  in  Salt  Lake  City  (Alder  et 
al.  1998). 

David  Hogan’s  examination  of  30  years’  worth  of  Salt  Lake  Valley  dense 
fog  events  identified  the  critical  role  that  inversions  played  in  the  formation  of 
these  events.  This  study  identified  two  key  criteria  that  an  inversion  should 
exhibit  to  allow  dense  fog  formation;  the  requirement  that  an  inversion  be 
relatively  shallow,  less  than  3,700  feet  in  height,  and  an  inversion  have  moderate 
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strength,  exhibiting  greater  than  1.5°C/1000  feet  lapse  rate.  Without  meeting 
these  criteria,  the  formation  of  dense  fog  at  KSLC  is  very  unlikely  (Hogan  1998). 

In  2004,  Jonathan  Slemmer,  a  meteorologist  with  the  National  Weather 
Service  (NWS)  Center  Weather  Service  Unit  (CWSU)  in  Salt  Lake  City, 
completed  an  unpublished  study  on  dense  fog  at  KSLC.  The  climatological 
database  used  for  the  study  was  the  surface  observations  taken  at  KSLC  from 
July  1971  through  June  2001.  Slemmer  sought  to  investigate  the  frequency  and 
intensity  of  dense  fog  events  at  the  airport  and  what  the  impact  of  the  events  had 
on  aviation.  For  this  study  Slemmer  defined  “dense  fog”  as  “reported  surface 
visibility  equal  to  or  less  than  one-quarter  statute  mile  (SM)  without  precipitation 
other  than  snow  grains”;  a  definition  which  will  be  maintained  throughout  this 
paper. 


In  the  study,  the  author  clearly  defines  the  meteorological  parameters  that 
are  typically  observed  during  dense  fog  events.  Three  of  the  most  unique  and 
possibly  consequential  findings  presented  by  Slemmer  indicate  that: 

•  Dense  fog  is  primarily  an  early-winter  event,  with  77%  of  all 
occurrences  observed  in  December  and  January.  February 
accounts  for  18%  of  all  events  while  only  2%  of  events  occurred  in 
November  or  March,  with  no  other  months  providing  statistically 
significant  occurrences.  While  the  author  does  not  speculate  on 
the  implication  of  this  information,  it  indicates  a  seasonally-driven 
forcing  mechanism  that  may  be  related  to  the  interaction  between 
colder  winter  air  masses,  increased  frequency  of  frontal  systems 
and  a  slowly  cooling  lake.  It  is  likely  that  the  moisture  contribution 
to  the  near  surface  levels  around  the  GSL  decreases  significantly 
as  the  lake  temperature  falls  through  the  winter.  In  addition  to  the 
greater  likelihood  of  ice  formation  as  the  temperatures  lowers, 
colder  water  has  lower  vapor  pressure  and  contributes  less  water  to 
the  atmosphere,  lowering  the  likelihood  of  dense  fog  formation  as 
the  winter  goes  on. 

•  Winds  from  the  west  through  northeast  (260°-  040°  in  degrees  true) 
and  southeast  (130°-190°  in  degrees  true)  provide  two  distinct 
maxima  for  the  observed  wind  direction  during  dense  fog  events 
(see  Figure  7).  While  the  NE-oriented  maximum  is  to  be  expected 
as  the  GSL  can  be  found  across  this  azimuth  range  and  winds  from 
this  direction  are  expected  to  be  moister,  the  SE-oriented  maximum 
may  be  related  to  cold  moist  air  draining  down  the  Jordan  River 
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Valley  from  Utah  Lake  to  the  southeast.  Another  possible 
explanation  for  the  southeasterly  maxima,  and  one  believed  to  be 
occurring  during  the  event  that  is  the  focus  of  this  study,  is  that  the 
southeasterly  flow  is  related  to  warm  air  advection  and 
strengthening  or  creation  of  an  inversion  over  KSLC  and  the 
persistent  cold  dome  over  the  GSL. 


Figure  7.  Dense  fog  occurrences  at  Salt  Lake  City  International  Airport  by  wind 

direction.  (From  Slemmer2004) 

•  Less  than  6%  of  all  dense  fog  cases  occurred  with  the  temperature 
above  freezing.  This  is  important  to  note  because  during  the  24 
January,  2009  fog  event  at  KSLC  the  temperature  remained  above 
freezing  for  the  duration.  This  implies  that  this  case  is  atypical  in 
this  regard  and  should  be  considered  in  the  future  if  other  dense 
fog  events  are  to  be  evaluated  through  ESA. 

In  an  examination  of  dense  fog  events  that  occurred  in  December  and 
January  during  the  winter  of  2001  and  2002,  Slemmer  categorized  the  dense  fog 
events  into  three  distinct  event  types:  a  prolonged  inversion  case;  precipitation 
falling  into  a  weakened  inversion  and  shallow  cold  pool  cases;  and  shallow  cold 
pool  advecting  from  GSL  cases.  It  is  believed  that  the  event  that  is  the  focus  of 
this  study  is  defined  by  the  both  the  second  and  third  event  type. 
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One  additional  phenomenon  described  by  the  author  is  how  shallow  cold 
pools  over  the  GSL  may  “slosh”  back  and  forth  in  the  valley,  leading  to  possible 
repeating  progression  and  retreat  of  cold  air  masses  over  KSLC.  Though  it  is  not 
examined  in  this  study,  because  of  the  intensity  and  duration  of  the  dense  fog 
event  of  24  January,  2009,  it  is  not  believed  that  this  type  of  phenomenon  is 
impacting  the  mesoscale  circulation  of  this  study. 

F.  ATMOSPHERIC  CHARACTERIZATION 

1.  Synoptic  Analysis 

The  week  leading  up  to  the  dense  fog  event  at  KSLC  was  characterized 
by  an  “omega”  blocking  pattern  with  strong  ridging  over  the  western  U.S.  and  an 
upper-level  ridge  axis  oriented  meridionally  over  the  Western  Rockies.  In  the 
Salt  Lake  City  area,  this  meant  clear  cold  nights  with  lows  near  20  °F  and  sunny 
days  with  highs  in  the  mid-30s°F.  Visibility  at  KSLC  remained  relatively  favorable 
during  this  period  with  haze  and  mist  reducing  visibility  to  between  3  SM  and 
5  SM  at  nearly  all  hours,  as  might  be  expected  with  a  stagnant  air  mass  in  an 
intermountain  lake  basin  such  as  the  Great  Salt  Lake.  One  period  of  significantly 
reduced  visibility  did  occur  on  the  night  of  20  January,  as  visibility  at  KSLC 
fluctuated  from  1/4  SM  to  3  SM  in  freezing  fog  and  mist,  as  a  shallow  fog  bank 
vacillated  over  the  airfield.  The  shallow  and  localized  nature  of  this  feature  can 
be  implied  through  observations  at  Hill  Air  Force  Base  (ICAO:  KHIF)  and  Ogden- 
Hinckley  Field  (ICAO:  KOGD).  These  locations  situated  25-30  SM  north  of  KSLC 
and  200-500  feet  higher  in  elevation,  remained  unrestricted  in  visibility 
throughout  the  same  time  period,  though  both  are  nearly  the  same  distance  from 
the  boundary  of  the  GSL.  This  event  is  interesting  to  note  because  it  serves  as 
the  only  time  period,  excepting  the  later  rain  and  dense  fog  event  which  is  the 
focus  of  this  study,  when  visibility  was  significantly  reduced  at  KSLC  during  the 
seven  days  that  the  western  U.S.  was  under  the  influence  of  this  regime. 

Beginning  on  20  January  and  continuing  through  23  January,  the  strong 
ridging  across  the  Rockies  began  to  slowly  weaken  as  the  upper-level  wave 
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pattern  over  the  U.S.  transitioned  from  highly  meridional  to  weakly  meridional. 
As  this  transition  occurred,  the  ridge  over  western  U.S.  maintained  its  north-south 
orientation  while  slowly  decreasing  in  amplitude  (see  Figure  8). 


NWVESR.  Phy$koi  5cWne»»  MvWon 


500mb  OeopofonGal  Height  (m)  Composite  Moon 
Dl/20/09 

NCEP/NCAR  Reoiwlyjis 


15800 


15700 


5600 


5500 


: 


5300 


I  5100 


5000 


Figure  8.  500  hPa  Geopotential  Height  (m)  Composite  Mean  change  from  20 

January  2009  to  23  January  2009.  (From  NOAA/ESRL  Physical 
Sciences  Division:  http://www.esrl.noaa.gov/psd/) 

Further  evidence  of  the  weakening  ridge  can  be  found  in  the  METAR 
observations  across  northern  Utah  from  21-23  January  which  indicate  increasing 
upper  and  mid-level  cloud  cover  and  lowering  altimeters.  Soundings  launched 
from  KSLC  at  this  time  show  decreasing  static  stability  aloft  and  increasing 
moisture  at  nearly  all  vertical  levels,  though  a  strong  low-level  inversion  remains 
present  throughout. 


By  23  January,  the  poleward  extent  of  the  ridge,  and  its  influence  over  the 

western-central  U.S.,  has  weakened  sufficiently  to  allow  a  mid-level  system  to 

push  across  the  Northern  Rockies.  As  this  system  moved  eastward,  it  was  led 

by  a  surface  cold  front  that  pressed  down  the  east  side  of  the  Rockies  and  into 

the  Great  Plains  early  on  23  January.  Behind  this  cold  front  a  strong  surface 

anticyclone  amplified  into  the  Great  Plains  during  the  day  on  23  January,  with 
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1047  hPa  central  pressure  estimated  over  southern  Alberta,  and  surface 
pressure  values  of  1040  hPa  reaching  into  northwest  South  Dakota  (see  Figure 
9).  At  the  same  time  that  the  cold  front  progressed  southward  across  the 
northern  Great  Plains,  a  surface  low  with  attendant  cold  front  entered  the 
southwest  U.S.,  pushing  into  southern  California.  This  disturbance  was  relatively 
weak  (lowest  central  pressure  recorded  as  1014  hPa  at  1200  UTC  on  23 
January)  and  its  eastward  progression  was  stopped  by  ridging  still  present  over 
the  Rockies. 


Figure  9.  Hydrometeorological  Prediction  Center  Surface  Analysis 

from  2100  UTC  on  23  January,  2009. 

(From  http://www.hpc.ncep.noaa.gov/) 


The  combined  effect  of  these  two  features  was  to  force  a  warm,  southerly 
flow  northward  across  the  southern  and  central  Intermountain  West  during  the 


31 


afternoon  hours  of  23  January,  2009.  Note  in  Figure  9,  the  southerly  winds 
recorded  at  stations  in  Arizona,  Nevada,  and  Utah,  and  the  generally  warm  air  in 
the  region  with  temperature  readings  of  73 °F  in  Phoenix,  61  °F  in  Las  Vegas,  and 
50 °F  in  southwestern  Nevada.  This  push  of  warm  air  over  the  cold  continental 
air  in  place  in  northern  Utah  was  essential  in  developing  conditions  ideal  for  the 
dense  fog  event  at  KSLC  early  on  24  January. 

2.  Precursor  Rain  Event 

On  23  January,  the  atmosphere  in  northern  Utah  had  become  sufficiently 
destabilized  with  the  passage  of  the  mid-level  disturbance  (Figure  10)  to  support 
vertical  motion  yielding  precipitation  in  the  form  of  light  rain  at  KSLC,  KHIF, 
KOGD,  and  as  far  south  as  Provo  (ICAO:  KPVU). 
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Figure  1 0.  Hydrometeorological  Prediction  Center  500  hPa  upper-air  analysis 
containing  observations,  heights,  and  temperatures  for  1200  UTC  23 
January  2009.  The  black  box  indicates  the  positioning  of  the  weak 
disturbance  moving  through  northern  Utah  on  23  January  2009. 

On  23  January,  Salt  Lake  City  Airport  recorded  0.11”  of  rain  from  0600 

UTC  to  1200  UTC  and  later  recorded  0.10”  of  rain  between  1600  UTC  and  2000 

UTC.  During  the  initial  rainfall  period,  surface  winds  at  all  of  the  local  stations 

were  generally  from  the  north,  though  winds  at  KSLC  exhibited  significant 

variability,  with  wind  source  alternating  90-180  degrees  every  other  hour,  while 

remaining  below  lOkts.  Winds  during  the  second  rain  event  are  observed  to  be 

predominantly  from  the  south,  though  significant  variability  is  still  observed.  Also 

important  to  note  is  the  steady  increase  in  temperature  on  23  January  at  KSLC, 

as  temperatures  rose  from  01  °C  at  1500  UTC  to  08  °C  at  2300  UTC,  while  light 

rain  fell  frequently.  Though  some  variability  in  wind  direction  should  be  expected 

in  a  geographically  complex  area,  the  significant  rise  in  temperatures  at  KSLC 
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and  high  variability  in  wind  direction  suggest  that  a  warm  frontal  boundary  had 
formed  in  the  southeast  GSL  basin.  This  boundary  developed  as  warm  air  from 
southern  Utah  met  the  colder  air  mass  over  the  GSL  Basin  and  the  influence  of 
the  mid-level  disturbance.  Where  these  two  air  masses  met  is  where  vertical 
motion  was  initiated  as  the  warm  air  rose  isentropically  up  and  over  the  colder 
dense  lake  air  (see  Figure  1 1 ). 


Figure  1 1 .  Radar  imagery  over  Salt  Lake  Basin  from  1 800  UTC 
23  January  2009.  (From  National  Climatic  Data  Center 
http  ://www.  ncdc.  noaa.gov/radar-data) 

The  variability  in  the  surface  observations  suggest  that  these  air  masses 

converged  over  the  southern  GSL  basin  for  several  hours  midday 

on  23  January,  before  the  warm  front  finally  displaced  the  cold  dome  over  KSLC 

between  1800  UTC  and  1900  UTC.  After  the  warm  front  moved  past  KSLC 

temperatures  rose  quickly,  precipitation  ended,  and  visibility  increased  to  3  SM. 

The  1200  UTC  23  January  KSLC  skew-T  (Figure  12)  indicates  a  strong  low-level 
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inversion  still  present,  with  southerly  winds  above  the  inversion,  veering  to 
westerly  at  750  hPa.  Examination  of  the  0000  UTC  24  January  skew-T  (Figure 
13)  shows  that  enough  destabilization  had  occurred  with  the  passage  of  the  mid¬ 
level  disturbance  to  mix  out  the  low-level  inversion  that  had  been  in  place  since 
17  January.  Note  that  the  low  level  winds  at  0000  UTC  have  become  northerly, 
which  coincides  with  the  onset  of  the  dense  fog  event. 


Figure  1 2.  1 200  UTC  23  January  2009  sounding  from  KSLC  indicating 

presence  of  strong  low-level  inversion. (From  University  of  Wyoming, 
http://weather.uwyo.edu/upperair/sounding.html) 
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Figure  1 3.  0000  UTC  24  January  2009  sounding  from  KSLC  indicating  that  the 
low-level  inversion  has  been  ameliorated.  (From  University  of 
Wyoming,  http://weather.uwyo.edu/upperair/sounding.html) 

By  1200  UTC  on  24  January  (Figure  14),  the  low-level  inversion  has 
reformed  due  primarily  to  a  decrease  in  temperatures  at  the  near-surface  level, 
as  temperatures  at  the  top  of  the  boundary  layer  remained  unchanged.  With  the 
inversion  reformed,  southerly  and  southeasterly  flow  returned  as  the  dense  fog 
event  concluded. 
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Figure  1 4.  1 200  UTC  24  January  2009  sounding  from  KSLC  indicating  the 
reformation  of  the  low-level  inversion. (From  University  of  Wyoming, 
http://weather.uwyo.edu/upperair/sounding.html) 

Though  no  skew-T’s  are  available  directly  over  the  lake,  it  is  hypothesized 

that  that  an  inversion  remained  in  place  through  the  duration  of  this  event  with 

surface  water  temperatures  across  the  lake  in  the  30s  °F  (Adams),  and  due  to  the 

relative  weakness  of  the  mid-level  disturbance.  Though  the  low-level  inversion  at 

KSLC  mixed  out  by  0000  UTC  on  24  January,  the  strength  of  the  inversion  over 

the  lake  was  likely  enhanced  after  the  warm  southerly  winds  met  the  boundary 

with  low-level  cold  pool  over  the  GSL.  The  isentropic  uplift  over  the  cold  dome 

that  ensued  strengthened  the  inversion  over  the  lake  and  forced  air  parcels  up 

and  over  the  inversion,  inducing  precipitation  which  fell  through  the  inversion 

toward  the  lake  surface,  fully  saturating  an  already  moist,  cool  environment. 

Adding  to  the  low  level  moisture  across  the  lake  was  widespread  rainfall  that  fell 

across  the  GSL  earlier  in  the  day  during  the  passage  of  the  mid-level  disturbance 
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(not  shown).  This  type  of  fog  formation  was  identified  by  Hogan  as  a  rare  but 
definable  case  when  light  precipitation  falling  into  an  already  moist  stable  layer 
causes  the  formation  of  dense  fog  (Hogan  1988).  As  the  warm  air  advection 
stopped  with  the  weakening  low  to  the  southwest,  the  cold,  moist,  and  dense  air 
over  the  GSL  flows  southward  enveloping  the  southern  GSL  basin  in  dense  fog. 

3.  Dense  Fog  Event 

At  2243  UTC  on  23  January,  winds  at  KSLC  were  calm,  visibility  was  3 
SM  in  haze,  and  the  temperature  and  dew  point  were  08C  and  05C  respectively. 
Within  one  hour,  conditions  changed  considerably  as  winds  became  northerly  at 
06  kts  with  visibility  lowering  to  1/16  SM  in  dense  fog  and  temperature  falling  to 
04 °C  at  2342  UTC.  Sunset  locally  on  23  January,  2009  was  at  0033  UTC  (24 
January)  or  1733  MST  (23  January),  and  occurred  1-2  hours  after  the  onset  of 
the  dense  fog  event.  Though  diurnal  effects  may  have  partially  contributed  to  the 
development  of  this  event,  cloud  cover  over  KSLC  remained  broken  to  overcast 
throughout  the  day  on  23  January  limiting  insolation,  so  it  is  unlikely  that  the 
daytime  heating  played  a  substantial  role  in  the  formation  of  this  event.  The 
average  high  temperature  at  KSLC  for  23  January  is  38  °F,  while  the  temperature 
recorded  this  date  was  48  °F,  indicating  substantial  anomalous  warming.  Dense 
fog  (<  1/16  SM  or  less)  and  light  (<06  kts)  west-to-northwesterly  or  calm  winds 
remained  predominant  at  KSLC  for  6  hours  following  onset  of  the  event,  through 
05  UTC  on  24  January.  The  temperature  at  KSLC  dropped  to  01  °C  at  0153  UTC 
on  24  January.  At  0553  UTC  on  24  January  (2253  MST,  23  January),  the  first 
sign  of  a  break  in  the  event  is  observed  as  winds  shifted  to  southeast  while 
remaining  light  (<  5  kts)  and  visibility  improved  slightly  to  1/8  SM.  These 
conditions  prevailed  until  1044  UTC  when  the  southeasterly  flow  strengthened  to 
7  kts  and  visibility  improved  to  3  SM.  Within  a  few  minutes  (1053  UTC)  the 
southeasterly  winds  were  measured  at  8  kts  and  visibility  rose  to  5  SM  as  the 
dense  fog  event  rapidly  eroded.  The  total  time  spent  with  visibility  at  or  near  1/8 
SM  is  approximately  10  hours. 
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4.  Relation  of  Physical  Interpretation  to  ESA 

The  importance  of  the  cold  dome  over  the  GSL  in  the  creation  and 
maintenance  of  low  level  inversions  is  believed  to  be  a  primary  factor  in  the 
development  of  fog  at  Salt  Lake  City  International  Airport  for  this  event.  It  is 
following  the  breakup  and  subsequent  reforming  of  the  inversion  that  the  cold 
moist  air  over  the  lake  is  able  to  flow  southward  over  KSLC.  It  is  believed  that 
the  breakup  of  the  inversion  is  enabled  by  the  vertical  forcing  associated  with  of 
the  mid-level  trough  and  the  advancing  warm  front. 

The  mesoscale  system  that  develops  within  the  synoptic  forcing  provides 
an  initial  predictor  as  to  which  meteorological  variables,  and  which  locations,  may 
be  important  in  the  formation  of  dense  fog  at  KSLC.  The  physical  interpretation 
described  in  the  atmospheric  characterization  identify  the  potential  role  that  the 
surface  and  planetary  boundary  layer  (PBL)  temperature,  wind,  and  relative 
humidity  may  offer  as  useful  variables  in  the  development  of  a  complete  ESA  for 
this  region. 

G.  PARAMETERS 

Previous  work  by  Chilcoat  (2012)  identified  the  need  for  careful  selection 
of  both  analysis  state  (x)  and  forecast  ( J )  parameter  to  determine  the  most 
responsive  and  indicative  sensitivities.  Chilcoat  found  the  u  and  v-wind 
components  on  the  lowest  model  layer  (~20  meters)  to  be  the  most  reliable 
analysis  state  parameter  for  predicting  robust  sensitivities.  For  forecast 
parameter,  Chilcoat  averaged  ensemble  average  water  vapor  mixing  ratio,  qv,  in 
a  2x2  (8-km  x  8-km)  square,  hereafter  referred  to  as  a  J-box,  located  over  KSLC 
and  summed  over  the  first  two  model  levels  (~58  meters). 

1.  Type  Selected  For  Study 

This  study  instead  found  the  most  robust  sensitivities  with  temperature  as 
the  analysis  state  parameter.  While  u  and  v  component  winds  did  exhibit 
sensitivity  in  our  forecast  system,  the  magnitude  was  -25%  of  that  found  with 
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temperature.  Other  state  parameters  were  also  considered,  including  qv  and 
total  dry  air  mass  in  the  column  (p).  Neither  parameter  provided  a  reliable, 
robust,  or  coherent  sensitivity  relationship  in  our  forecast  system. 

2.  Configuration  of  Parameters 

Alterations  were  also  necessary  to  the  size  of  the  J-box  centered  over  the 
airport,  though  qv  was  remained  the  forecast  parameter.  Due  to  the 
implementation  of  the  linearity  filtering  described  in  section  3.B.,  the  2x2  J-box 
summed  over  the  lowest  two  model  levels  used  by  Chilcoat  no  longer  provided  a 
reliable  sensitivity  response.  Sensitivities  using  Chilcoat’s  configuration 
continued  to  provide  a  reliable  sensitivity  response  at  the  6-hour  forecast  time, 
however,  12-hour  forecasts  using  this  configuration  failed  the  Gaussianity  testing. 
This  failure  was  determined  by  Lilliefors  test  results  indicating  that  the  ensemble 
distribution  of  forecasted  qv  values  in  the  area  enclosed  by  the  J-box  was  non- 
Gaussian. 

In  order  to  ensure  Gaussianity  of  forecast  parameter  and  continuity 
between  forecast  lead  times,  it  was  necessary  to  alter  the  size  and  position  of  the 
J-box.  Multiple  combinations  of  J-box  sizing  were  attempted,  including  reduction 
of  the  J-box  to  a  single  point  closest  to  KSLC,  as  well  as  expansion  to  horizontal 
grid  spacing  of  3x3  and  4x4,  with  variations  of  model  levels  between  level  1 ,  level 
2,  and  the  summation  of  these  levels.  It  was  also  theorized  that  expansion  of  the 
J-box  toward  the  expected  source  (the  GSL)  would  provide  a  more  reliable 
response.  Due  to  the  nearly  equal  proximity  of  the  GSL  to  the  north  and  west, 
and  also  the  variation  of  the  low  level  winds  leading  up  to  the  dense  fog  event,  it 
was  speculated  that  expansion  of  the  J-box  outward  in  all  directions  would 
provide  the  most  representative  data  field  for  low  level  q^  around  the  airport.  . 
After  trial  and  error,  the  smallest,  yet  still  Gaussian,  J-box  that  could  be 
determined  to  be  valid  at  the  6-hour  and  12-hour  forecast  was  a  3x3  box  (12  km 
x  12  km)  centered  nearly  directly  over  KSLC,  and  summed  over  the  first  two 
model  levels.  The  selection  of  a  larger  J-box  than  was  previously  used  by 
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Chilcoat  is  should  improve  ESA  findings  as  the  area-averaging  of  J  provides 
more  robust  results  by  averaging  out  some  of  the  sample  error  and  noise 
associated  with  the  distributions  (G.  J.  Hakim  2013,  personal  communication). 

Though  the  specific  dimensions  of  Chilcoat’s  J-box  were  not  retained  for 
this  study,  this  does  not  invalidate  the  results  of  the  previous  work.  Linearity 
testing  conducted  by  Chilcoat  (2012,  section  4.D.)  determined  that  the  sensitivity 
results  he  achieved  were  linear.  This  allows  the  conclusion  to  be  drawn  that  the 
ensemble  distribution  of  x  and  J  he  selected  were  approximately  Gaussian.  The 
Gaussianity  testing  in  this  study  provides  an  a  priori  means  of  ensuring  that 
observed  and  forecast  parameters  are  drawn  from  normal  distributions,  which 
enables  implementation  of  ESA  with  greater  certainty. 

To  ensure  that  x  was  also  drawn  from  a  normal  distribution,  Lilliefors 
testing  was  also  applied  to  the  ensemble  distribution  of  temperature  at  the  most 
sensitive  grid  point  (Y  =  52,  X  =  121)  in  the  innermost  domain.  These  results  (not 
shown)  indicated  that  analysis  values  for  the  1800  UTC  23  January  2009 
ensemble  system  (6-hour  forecast)  were  Gaussian.  The  results  of  testing  for  the 
1200  UTC  23  January  ensemble  system  (12-hour  forecast)  were  mixed.  While 
Gaussianity  tests  at  the  most  sensitive  point  were  weakly  non-Gaussian,  tests  at 
several  adjacent  sensitive  grid  points  (Y=51,  X=121  and  Y=53,  X=120)  with 
similar  levels  of  sensitivity  indicated  entirely  Gaussian  distributions  of 
temperature  (see  figure  15).  This  finding  implies  that  ensemble  distribution  of 
temperature  at  the  concentrated  area  of  greatest  sensitivity  is  effectively  normal. 
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a.  b. 


Figure  1 5.  Histograms  of  ensemble  analysis  potential  temperature  (6)  for  the 
lowest  model  level  at  1200  UTC  23  January  2009  for  (a.)  the  most 
sensitive  grid  point  in  the  innermost  domain,  (Y=52,  X=121)  and  (b.) 

an  adjacent  grid  point  (Y=51 ,  X=1 21 ).  The  red  line  indicates  the 
normal  distribution  while  the  blue  bars  indicate  the  probability  density 
function  of  6  at  each  grid  point. 

The  conclusion  that  can  be  drawn  from  these  tests  is  that  ensemble 
distribution  of  the  forecast  parameter  (qv)  at  KSLC  is  Gaussian,  while  the 
analysis  state  variable  (7)  at  the  most  sensitive  grid  point,  is  Gaussian,  or  that  it 
represents  the  normal  distribution  with  sufficient  closeness  to  meet  the  required 
assumptions  necessary  to  perform  the  linear  regression  that  forms  the  basis  of 
ESA. 
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IV.  RESULTS 


A.  WRF  ENSEMBLE  PREDICTIONS 

Though  the  purpose  of  this  ensemble  study  was  not  to  quantify  accuracy 
of  the  WRF  mesoscale  ensemble  system,  the  qualitative  performance  of  the 
system  is  of  great  relevance.  To  accurately  identify  the  statistical  and  dynamic 
relationships  of  this  system,  it  is  necessary  for  the  ensemble  system  used  to 
define  the  system  to  have  skill.  This  section  will  discuss  the  performance  of  the 
96-member  WRF-ARW  model  used  for  this  experiment. 

1.  Moisture  in  GSL  Basin 

Accuracy  of  the  forecast  system  was  estimated  using  predictions  for  a  grid 
point  representing  KSLC  (X=137,  Y=88).  This  grid  point  lies  approximately  0.6 
miles  (1 .0  kilometer)  northeast  of  the  airport  and  is  most  representative  of 
conditions  at  KSLC. 

Figure  16  shows  EM  qv  (in  kg  kg"1)  at  a]_i  (the  lowest  model  level)  in  3-hour 
intervals,  clockwise  from  top  left,  from  1800  UTC  23  January  to  0600  UTC  24 
January.  This  figure  represents  the  96-member  ensemble  initialized  at  1800 
UTC  23  January.  The  image  shows  the  strong  concentration  of  water  vapor  in 
the  lowest  model  level,  represented  by  warmer  colors,  over  KSLC  at  1800  UTC 
23  January.  It  appears  also  at  this  time  that  stronger  concentrations  of  water 
vapor  extend  beyond  the  boundary  of  the  GSL,  nearly  overcoming  KSLC.  At 
2100  UTC,  the  higher  water  vapor  concentrations  have  expanded  further  from 
GSL,  clearly  encompassing  KSLC.  The  densest  water  vapor  concentrations 
push  outward  from  the  GSL,  filling  the  valleys  south  of  lake  10-20  kilometers 
inland.  Expansion  of  the  higher  water  vapor  concentrations  continues  through 
0000  UTC  on  24  January. 

Between  0000  UTC  and  0300  UTC  24  January,  the  horizontal  extent  of 

the  strongest  concentrations  of  water  vapor  begins  to  recede,  indicating  drying 

across  GSL  basin.  The  amelioration  of  the  dense  water  vapor  values  intensifies 
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between  0300  UTC  and  0600  UTC  24  January.  It  appears  that  KSLC  is 
surrounded  by  a  significant  dry  slot  by  0600  UTC  24  January,  indicating  that  the 
dense  fog  event  has  ended. 


Figure  16.  Ensemble  mean  (in  kg  kg"1)  at  r/-i  in  3-hour  intervals,  clockwise 
from  top  left,  from  1800  UTC  23  January  to  0600  UTC  24  January. 
UTC  annotated  as  Z.  96-member  ensemble  initialized  at  1800  UTC 
23  January.  Approximate  position  of  KSLC  is  represented  by  black 
star.  Warm  colors  represent  stronger  concentration  of  water  vapor. 

2.  Moisture  at  KSLC 

Moisture  at  KSLC  is  measured  using  an  Automated  Surface  Observing 
System  (ASOS)  at  a  height  of  approximately  2  meters.  To  compare  the  actual 
observed  amount  of  water  vapor  at  KSLC  versus  the  model  forecast  value, 
calculations  were  made  to  convert  the  observed  measurements  of  temperature, 
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dew  point  and  pressure  into  water  vapor  mixing  ratio,  in  units  of  kg  kg"1.  Water 
vapor  mixing  ratio  is  one  of  the  five  primary  output  fields  produced  by  WRF-ARW 
used  to  track  the  different  state  parameters  of  water. 

Figure  17  shows  the  forecasted  20-meter  and  2-meter  q„  at  KSLC  versus 
observed  at  KSLC.  The  forecast  values  represent  the  12-hour  ensemble 
mean  (EM)  forecast  values  from  the  ensemble  that  was  initialized  at  1800  UTC 
23  January.  Forecasts  from  this  system  end  at  0600  UTC  on  24  January  2009. 
Both  20-meter  and  2-meter  forecast  were  plotted  to  show  low-level  moisture 
continuity  at  near-surface  levels.  The  figure  indicates  that  at  1800  UTC  23 
January,  the  ensemble  system  is  already  nearing  its  12-hour  maximum  in  water 
vapor  mixing  ratio  concentration.  Both  the  2-meter  and  20-meter  water  vapor 
levels  increase  until  2100  UTC  23  January,  when  both  reach  maximum  levels. 
These  maximum  values  are  maintained  until  0300  UTC  on  24  January  when  the 
water  vapor  levels  begin  to  decrease. 
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Figure  1 7.  Ensemble  mean  forecast  20-meter  and  2-meter  at  KSLC  versus 
observed  qKat  KSLC.  The  red  solid  line  represents  20-meter  q^for 
the  forecast  initialized  at  1 800  UTC  23  January  while  blue  solid  line 
represents  2-meter  qv  from  the  same  initialization.  The  q^  recorded 
at  the  airport  is  represented  by  the  blue  dashed  line.  Both  forecasts 
end  at  0600  UTC  24  January. 


The  values  presented  in  Figure  17  show  a  marked  difference  in  both 
magnitude  and  timing  between  the  ensemble  forecast  and  the  observed  values. 
Though  it  would  be  preferable  to  have  the  magnitude  of  the  forecast  values  more 
closely  represent  the  actual  observations,  the  general  agreement  between  the 
two  (both  indicating  high  levels  of  water  vapor)  is  encouraging.  Observed  water 
vapor  mixing  ratios  during  the  event  fluctuated  significantly  from  3.91x1  O'3  kg  kg"1 
to  4.74x1  O'3  kg  kg'1.  Water  vapor  remained  high  even  as  the  dense  fog  lifted. 
Interestingly,  water  vapor  mixing  ratios  remained  high  at  4.42x1  O'3  kg  kg"1  at 
1200  UTC  24  January,  even  as  surface  visibility  had  increased  to  7  SM.  In 
addition,  light  some  rain  was  reported  as  the  fog  lifted. 

More  important  than  the  magnitude  difference  in  the  water  vapor  mixing 
ratio  is  the  timing  of  the  increases  and  decreases.  The  trend  represented  in 
Figure  17  shows  a  notable  decrease  in  forecast  value  of  qu  between  0300  UTC 
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and  0600  UTC  24  January.  This  can  be  interpreted  as  the  end  of  the  dense  fog 
event  at  KSLC.  Observations  show  that  the  dense  fog  did  not  lift  until  0900  UTC 
24  January,  which  means  that  the  forecast  is  3-6  hours  early  in  its  dissipation. 

Though  forecast  concentration  of  qv,  and  exact  timing  of  its  dissipation  are 
not  entirely  well-represented,  Figures  16  and  17  do  show  agreement  in  meso- 
and  local-scale  of  the  process  of  early  dissipation  of  the  dense  fog.  Chilcoat 
(2012)  found  similar  results  in  his  study  and  his  conclusion  that  the  timing  of  the 
event  was  not  critical  to  the  effectiveness  of  the  ESA  will  be  echoed  here.  The 
important  characteristics  of  the  forecast  system  are  that  it  adequately  predicts  the 
presence  of  high  levels  of  moisture  across  the  GSL  Basin  and  at  KSLC,  and  that 
it  faithfully  represents  the  dissipation  of  the  fog  event,  both  of  which  are  evident 
in  this  study. 

B.  QUANTITATIVE  RESULTS 

ESA  was  conducted  using  the  methodologies  outlined  in  section  III.  The 
linearity  (Lilliefors)  testing,  ESA  calculation  and  statistical  sampling  error 
corrections  were  performed  in  order,  with  results  defined  below.  Temperature 
was  chosen  as  the  initial  state  variable  {xj)  and  water  vapor  mixing  ratio,  qv, 
selected  as  the  forecast  parameter  (J).  In  this  analysis,  temperature,  T  (K)  is 
determined  at  the  analysis  lowest  eta-level  (approximately  25  m  AGL)  while 
water  vapor  mixing  ratio,  qv  (kg  kg"1)  is  a  summation  of  the  lowest  two  model 
levels  (approximately  60  m  AGL)  in  the  area  of  the  12-km  x  12-km  J-box 
surrounding  KSLC  defined  in  section  III.G. 

Sensitivity  results  (c/J/c/x,)  quantify,  in  kg  kg'1  K"1,  how  changes  in  J  relate 
to  changes  in  x,.  Positive  dJ/dx,  means  that  J  increases  with  an  increase  in  xi; 
and  J  decreases  with  a  decrease  in  x,.  Negative  dJ/dxi  means  that  J  increases 
with  a  decrease  in  x,  and  vice  versa.  In  this  study  we  will  focus  on  positive 
change  in  dJ/dXj,  as  this  study  seeks  to  understand  how  state  variable 
temperature  relates  to  dense  fog  formation  at  a  later  time. 
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The  results  described  in  this  section  are  based  on  a  96-member  ensemble 
system,  though  experiments  involving  varying  sizes  of  ensemble  configurations 
were  performed  and  will  be  discussed  later.  Experiments  were  conducted  using 
both  6-hour  forecast  systems  (initialized  1800  UTC,  23  January  2009)  and  12- 
hour  forecast  systems  (initialized  1200  UTC,  23  January  2009).  White  space  in 
Figures  18  and  21  indicate  where  data  has  been  masked  at  the  95%  confidence 
interval.  The  star  next  to  J  indicates  the  approximate  location  of  the  Salt  Lake 
City  International  Airport. 

1.  ESA  (6  hour) 

Results  from  the  6-hour  ESA  are  displayed  in  Figure  18.  Region  1,  which 
lies  approximately  150km  south/southwest  of  the  GSL,  indicates  the  most 
positive  sensitivity  region  within  the  innermost  domain.  Sensitivity  in  this  area 
ranges  from  4.433x1 O'4  kg  kg"1  K"1  to  6.61 5x1  O'4  kg  kg'1  K"1.  At  the  most 
sensitive  grid  point  (39.5072N  112.9197W)  a  ±0.0516  K  (1  oT)  change  in 
temperature  at  the  lowest  model  level  will  lead  to  an  equivalent  ±1.01 5x1  O'4  kg 
kg"1  (1  oqv)  change  in  the  6-hour  qv  at  KSLC.  For  comparison,  it  is  noted  that  this 
value  approximates  the  2m  and  20m  ensemble  mean  forecast  variation  in  qv  over 
the  same  time  period  (see  Figure  17). 
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Figure  18.  Sensitivities  (dJ/dxj)  using  1800  UTC  23  January  analysis  potential 
temperature  (in  K)  as  xh  and  forecasted  qv  (in  kg  kg"1)  as  J,  valid  at 
0000  UTC  24  January.  Units  are  kg  kg'1  K"1 

Other  strong  positive  sensitivity  responses  appear  in  Figure  18  across  the 
Central  Utah  Desert  Ranges,  indicated  at  Region  2  and  Region  3.  The 
geographical  commonality  between  these  three  regions  is  that  they  each  are 
located  near  elevated  terrain  features  that  have,  at  least  partially,  an  east-to-west 
orientation  (Figure  19).  In  addition,  these  features  are  positioned  where  a  strong 
plume  of  warm  air  pushes  northward  through  the  Central  Utah  Desert  on  23 
January  (see  Figure  20). 


49 


Figure  1 9.  WRF  terrain  map  indicating  elevations  between  1 000m  and  2000m. 

The  black  circle  (#1)  identifies  the  elevated  terrain  near  the  grid  point 
where  the  strongest  6-hour  ESA  sensitivity  response  is  located. 
Other  numbers  (#2-7)  indicate  other  regions  important  to  the 
sensitivity  response. 


Figure  20.  Ensemble  mean  6  (K)  at  q.-\  at  1800  UTC  23  January.  The  most 
sensitive  6-hour  ESA  grid  point  is  marked  in  the  black  circle  (#1), 
near  the  strongest  concentration  of  warm  air,  south  of  the  GSL. 
Other  numbers  (#2-7)  indicate  other  regions  of  sensitivity 

importance. 
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These  ranges  stand  in  contrast  to  the  largely  north-to-south  orientation  of 
much  of  the  terrain  across  northern  Utah,  as  discussed  in  section  III.E.  The 
alignment  of  these  features  is  important  because,  in  contrast  to  terrain  oriented 
from  north-to-south,  these  topographic  features  may  block  or  divert  southerly 
flow.  As  southerly  flow  in  this  region  is  likely  to  be  correlated  with  warm-air 
advection,  these  east-to-west  terrain  features  are  the  locations  where  warm-air 
advection  across  the  domain  is  indicated. 

A  different  sensitivity  response  is  evident  in  Region  4,  northeast  of  the 
GSL.  This  region  indicates  a  low-lying  plain  between  the  GSL  and  the  Wasatch 
Range.  The  sensitivity  evident  here  is  created  as  southerly  warm-air  advection  is 
forced  between  the  higher  terrain  of  the  Wasatch  and  the  cold  dome  over  the 
GSL.  Warm  air  spills  into  this  plain  northeast  of  the  lake,  creating  a  layer  of 
positive  sensitivities. 

Several  areas  of  strong  negative  sensitivity  are  also  present  in  Figure  18, 
though  their  magnitude  is  less  than  the  positive  regions,  with  no  values 
exceeding  -2.5x1 0"4  kg  kg"1  K'1.  The  negative  regions  indicated  on  the  map, 
Region  5  (Tooele  Valley),  Region  6  (southeast  Salt  Lake  Valley),  and  region, 
Region  7  (lee  side  of  the  Uinta  Mountains),  represent  low  lying  areas  with 
blocked  flow  where  cold  dense  air  is  likely  to  pool  during  stable  regimes. 

Analyzed  together  these  regions  of  positive  and  negative  sensitivity 
reinforce  the  conceptual  model  defined  in  Section  III.  This  model  identifies  the 
significance  of  the  warm-air  advection  prior  to  the  event,  as  well  as  the 
importance  of  the  low  level  stability  profile  that  supported  the  genesis  of  the 
dense  fog  event. 

2.  ESA  (12  hour) 

The  results  from  the  12-hour  ESA  in  Figure  21  show  a  somewhat  less 
clearly  defined  sensitivity  relationship  than  the  6-hour  ESA.  This  ESA  used  initial 
condition  temperature  from  1200  UTC  23  January,  and  forecast  qvat  0000  UTC 
24  January. 
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Figure  21 .  Sensitivities  ( dJ/dxi)  using  1200  UTC  23  January  analysis  potential 
temperature  (in  K)  as  x-h  and  forecasted  qv  (in  kg  kg"1)  as  J,  valid  at 
0000  UTC  24  January.  Units  are  kg  kg'1  K'1. 

The  strongest  positive  sensitivity  is  located  once  again  in  the  Central  Utah 
Desert  Range,  as  identified  by  Region  1  in  Figure  21.  This  region  exhibits  a 
sensitivity  range  from  3.82x1  O'4  to  4.1 93x1  O'4  kg  kg'1  K'1,  which  is  weaker  in 
magnitude,  but  still  geographically  close  to  the  strong  positive  sensitivity  region 
evident  for  the  6-hour  ESA.  Analysis  at  the  same  grid  point  examined  for  the  6- 
hour  ESA  (39.5072N  112.9197W)  indicates  a  ±0.0684  K  (1aT)  change  in 
temperature  at  the  lowest  model  level  will  lead  to  ±1.787x1  O'4  kg  kg"1  change  in 
the  12-hour  qvat  KSLC.  We  note  that  this  value  is  an  increase  (+  0.772x1 0"4  kg 
kg'1)  in  predicted  qv  from  the  6-hour  ESA.  This  result  is  important  because  it 
shows  a  continued  strong  sensitivity  relationship  across  a  specific  12-hour 
interval  which  Chilcoat  (2012)  identified  as  a  limitation  to  his  ESA  because  of  the 
nonlinear  dynamics  that  were  introduced  by  a  wind  shift. 
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The  other  strong  sensitivity  regions  in  the  Central  Utah  Desert  Range 
previously  evident  in  the  6-hour  ESA  (Figure  19,  Regions  2  and  3)  fail  the 
confidence  interval  testing.  This  change  suggests  that  the  warm  air  advection 
predominant  in  this  region  during  the  6-hour  ESA  is  a  weaker  predictor  of  Jover 
the  12-hour  interval.  The  continued  presence  of  this  strong  positive  sensitivity 
region  identifies  it  as  the  most  important  region  of  sensitivity  in  the  inner  grid. 

Also  notable  in  the  12-hour  ESA  are  the  positive  sensitivities  at  Region  2, 
3,  and  4,  evident  in  Figure  21.  Region  2  is  centered  over  the  Stansbury 
Mountain  Range,  which  rises  over  11,000  feet.  This  positive  sensitivity  over  an 
elevated  terrain  region  indicates  an  increase  in  vertical  propagation  of  warm-air 
advection  from  the  6-hour  ESA.  This  difference  in  vertical  distribution  of  positive 
sensitivity  indicates  the  influence  of  warm  air  at  higher  levels,  likely  enabled  by 
the  destabilizing  effects  of  the  mid-level  disturbance  which  allows  for  the  warm 
air  to  rise  further  in  the  atmosphere.  Region  3  is  very  close  to  a  6-hour  ESA 
negative  sensitivity  point  in  the  south  Salt  Lake  Valley  (Figure  18  -  Region  6). 
This  region  exhibits  a  sensitivity  sign  change  from  6-hour  to  12-hour  forecast  on 
a  very  short  spatial  scale  (~20km  across  the  Salt  Lake  Valley)  and  emphasizes 
this  area  as  important  to  the  understanding  of  the  dynamic  profile  of  this  event. 
Region  3  is  unique  in  that  it  is  blocked  on  three  sides  by  mountains  (Oquirrh 
Mountains  to  the  west,  Traverse  Range  to  the  south,  and  Wasatch  Range  to  the 
east),  yet  it  also  appears  to  serve  as  a  primary  corridor  for  warm-air  advection. 
Region  4  represents  the  Green  River  Valley  and  warm  sensitivities  here  can  be 
attributed  to  the  destabilizing  influence  of  the  mid  level-disturbance  which  helped 
to  mix  out  the  inversion  and  the  cold  dense  air  trapped  underneath. 

Also  in  Figure  21  are  weak  negative  sensitivities  west  of  the  GSL  and  over 
the  south  end  of  the  GSL.  These  sensitivity  responses  may  indicate  the 
influence  of  the  cold  dense  air  of  the  lake  as  it  expands  outward  from  the  GSL 
during  the  dense  fog  event. 
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3. 


Linear  Correlation 


To  further  examine  the  relationship  between  x,  and  J,  scatterplots  were 
created  for  from  6-hour  and  12-hour  forecast  data  sets.  These  plots  allow  simple 
examination  of  the  data  with  the  potential  to  identify  linear  (or  non-linear) 
relationships,  clustering  of  variables,  outliers,  and  changes  in  spread  of  each 
variable  as  a  function  of  the  other  (Wilks  2006).  Also  computed  was  correlation 
coefficient  between  x,  and  J,  labeled  as  r  on  the  scatterplots.  This  coefficient 
provides  a  measure  of  the  strength  and  magnitude  of  the  linear  relationship 
between  two  variables.  The  ESA  computed  sensitivity  value  ( dJ/dx )  also 
appears  on  each  scatterplot. 

Figures  22  and  23  represent  the  6-hour  and  12-hour  96-member 
ensemble  scatterplots  for  analysis  potential  temperature  at  1800  UTC  (6-hour, 
Figure  22)  and  1200  UTC  (12-hour,  Figure  23)  23  January  at  the  most  sensitive 
grid  point  (39.5072N  112.9197W)  versus  the  forecast  for  qv  at  KSLC  at  0000 
UTC  24  January.  The  most  sensitive  grid  point  lies  approximately  170  km 
southwest  of  KSLC  and  was  identified  as  Point  1  in  Figure  18. 

The  scatterplot  of  the  6-hour  forecast  data  suggest  good  linear  agreement 
between  the  two  parameters  with  some  spread  evident  and  a  minimal  number 
outliers.  The  correlation  coefficient  associated  with  this  data  is  0.55  and 
indicates  a  moderate  linear  correlation  between  x,and  J. 


54 


„„x  10 

3.6  r 


u 

_l 

* 


3.4 


Ol 

3  2 

> 

O 

T3 

<U 


3- 


r  =0.5465 

dJdx  =  6.61534e-04 


•• 


•• 


* 

»  *♦  «  « 

.•  .**•*.'  *•*.••  .  ••  . 

'  •  * ; 


2  0 _ i _ i _ i _ i _ i _ i _ i _ I 

285.65  285.7  285.75  285.8  285.85  285.9  285.95  286  286.05 

_ Analysis  Potential  Temnflratnrfi  < IQ  Rt  39  5072N  112  9197W _ 

Figure  22.  Scatterplot  of  1 800  UTC  23  January  analysis  9  at  the  most  sensitive 
grid  point  and  forecasted  qvat  KSLC  at  0000  UTC  24  January. 


Figure  23  represents  the  12-hour  forecast  for  the  same  parameters  and 
indicates  a  significant  breakdown  in  the  linearity  present  in  the  6-hour  forecast. 
The  associated  /"-value  of  this  scatterplot  is  r  =  0.023.  Analysis  of  the  scatterplot 
shows  increased  spread  between  parameters,  suggesting  the  relationship 
between  these  variables  is  significantly  less  linear  than  the  6-hour  forecast.  A 
closer  inspection  indicates  a  small  grouping  of  members  predicting  a  higher 
value  of  qv  (3.1 5x1 0"4  kg  kg"1  to  3.45x1  O'4  kg  kg'1)  than  the  majority  of  the 
ensemble  field.  This  grouping  of  10  members  is  the  likely  cause  for  the  low  r- 
value  and  suggests  that  there  is  another  predictor  for  J  at  this  same  initial  time 
(1200  UTC  23  January).  This  alternate  predictor  could  be  9  at  another  vertical  at 
the  same  sensitivity  spot,  9  at  a  separate  grid  point  and/or  /7-level,  or  an  entirely 
different  x,  variable. 
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Figure  23.  Scatterplot  of  1 200  UTC  23  January  analysis  potential  temperature 
at  the  6-hour  ESA  most  sensitive  grid  point  and  forecasted  qv  at 
KSLC  at  0000  UTC  24  January. 


The  significant  breakdown  in  linearity  between  6-hour  and  12-hour 
scatterplots  indicates  the  limitation  of  a  single  fixed  location  to  predict  sensitivities 
over  longer  (>6-hour)  time  intervals.  Examination  of  the  sensitivity  plots  at  both 
6-hour  (Figure  18,  Point  1)  and  12-hour  times  (Figure  21,  Point  1)  indicates  that 
the  most  sensitive  grid  point  moved  9  kilometers  to  the  northeast  from  the  6-hour 
forecast  to  the  12-hour  forecast.  In  addition,  the  clustering  of  positive  sensitivity 
points  around  the  maximum  sensitivity  grid  point  moved  notably  to  the  northeast 
and  formed  in  a  more  dispersed  geographic  distribution.  This  indicates  that  the 
most  sensitive  grid  point  for  one  forecast  time  scale  is  not  likely  to  be  the  most 
sensitive  grid  point  at  a  different  time  scale. 

To  justify  placement  of  future  sensors  using  the  results  of  an  ESA,  it  will 
be  necessary  to  create  an  observing  network  with  accurate  sensitivity  beyond  the 
6-hour  time  frame.  However,  the  limited  forecast  improvement  at  short  time 
scales  for  one  sensor  location  forces  the  consideration  of  a  network  of  sensors  at 
multiple  locations.  To  investigate  this  scenario,  sensitivities  and  linearity  were 
also  calculated  at  the  most  sensitive  and  second  most  sensitive  grid  point,  from 
data  provided  by  the  12-hour  ESA. 


56 


Figure  24  shows  scatterplots  of  1200  UTC  23  January  analysis  potential 
temperature  at  the  most  sensitive  (a.)  and  second  most  sensitive  (b.)  grid  points 
as  calculated  during  the  12-hour  ESA.  Both  analysis  potential  temperature 
distributions  are  plotted  against  forecast  qv  at  KSLC  at  0000  UTC  24  January. 
These  results  indicate  a  significant  increase  in  linearity  when  compared  to  the 
distribution  of  data  at  the  6-hour  ESA-indicated  most  sensitive  grid  point.  Neither 
/"-value  in  Figure  24  indicates  strong  or  even  moderate  linearity,  though  an  order 
of  magnitude  increase  is  evident  in  both  correlation  coefficient  and  sensitivity 
value.  Also  visible  in  Figure  24  is  the  contribution  from  an  alternate  predictor  that 
was  evident  in  Figure  23.  The  continued  presence  of  this  grouping  of  high  qv 
members  at  several  alternate  points  supports  the  proposition  that  the  influence  of 
another  predictor  is  inducing  non-linear  responses  over  this  12-hour  interval. 
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Figure  24.  Scatterplot  of  1 200  UTC  23  January  analysis  potential  temperature 
at  (a.)  the  most  sensitive  grid  point  from  the  1 2-hour  ESA  and  (b.) 
the  second  most  sensitive  grid  point  from  the  1 2-hour  ESA.  Both  are 
plotted  against  forecast  qvat  KSLC  at  0000  UTC  24  January. 


Also  indicated  along  the  x-axis  in  Figure  24  are  the  distances  in  kilometers 
between  the  6-hour  ESA  most  sensitive  point  and  the  12-hour  ESA  first  two  most 
sensitive  points.  Note  that  the  three  most  sensitive  grid  points  exist  across  a  30- 
kilometer  distance  in  varying  terrain.  These  data  suggest  that  a  future  observing 
network  in  complex  terrain  may  need  to  include  the  placement  of  multiple 
sensors  across  a  relatively  broad  swath  of  terrain. 

4.  Support  for  Conceptual  Model 

Evidence  supporting  the  above  conclusions  is  present  in  the 
synoptic/mesoscale  analysis  provided  in  Section  III  which  identified  the  role  that 


58 


the  mid-ievel  disturbance  played  during  the  day  on  23  January.  A  graphical 
representation  is  provided  in  Figure  25. 


Mid-level  disturbance  passes  N.  Utah  mid-day  on  23  Jan, 
disrupting  stability  profile  &  saturating  low  levels  over  GSL 
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Figure  25.  Schematic  interpretation  of  effects  of  23  January  2009  mid-level 
disturbance  and  synoptically  forced  warm-air  advection  that  led  to  a 
dense  fog  event  at  Salt  Lake  City  International  Airport  on  24  January 

2009.  (After  Google  Maps) 


A  schematic  representing  the  12-hour  ESA  profile  (Figure  26)  indicates 

the  complexity  of  the  forcing  across  this  time  interval.  The  mid-level  disturbance 

that  moved  through  northern  Utah  exhibited  greater  influence  across  the  GSL 

basin,  with  lesser  influence  felt  in  central  Utah.  This  disturbance  resulted  in 

vertical  mixing  of  the  stable  boundary  layer  over  northern  Utah,  which  likely 

would  have  increased  temperatures  in  low-lying  valleys  (positive  sensitivities 

annotated  in  Regions  3  and  4,  Figure  21)  that  had  filled  with  cold  dense  air 

beneath  low-level  inversions.  While  the  mid-level  disturbance  moved  through 
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northern  Utah,  warm-air  advection  was  also  ongoing  across  the  region.  Warm- 
air  advection  that  occurred  throughout  the  12-hour  interval  period  appears  at  the 
most  positive  sensitive  region  well  south  of  GSL  (Region  1 ,  Figure  21 ),  and  in  the 
elevated  terrain  immediately  south  of  GSL  (Region  2,  Figure  21).  Figure  26 
shows  a  schematic  of  the  1 2-hour  ESA. 
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Figure  26.  Conceptual  sensitivity  field  using  1 200  UTC  23  January  2009 
analysis  temperature  as  the  analysis  state  variable  (x,)  and  qvat 
0000  UTC  24  January  2009  as  the  forecast  metric  (J)  for  the  box 
centered  over  KSLC.  (After  Google  Maps) 


By  1800  UTC  on  23  January,  the  influence  of  the  mid-level  disturbance 

had  weakened  enough  that  warm  air  advection  began  to  define  different  6-hour 

sensitivity  relationships.  Strong  positive  sensitivities  were  evident  near  elevated 

terrain  features  with  at  least  some  east-to-west  aligned  orographic  profile 

(Region  1,  2,  and  3,  Figure  18)  and  where  warm  air  advection  was  being 
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tunneled  (Region  4,  Figure  18).  This  profile,  combined  with  strong  negative 
sensitivities  in  low  lying  areas  (Region  5,  6  and  7,  Figure  18)  provide  evidence 
that  this  flow  was  rebuilding  stability  into  northern  Utah.  The  resulting  inversion 
helped  support  the  advection  of  cold  moist  air  over  KSLC  which  defined  the 
dense  fog  event.  A  schematic  representing  these  6-hour  sensitivities  appears  in 
Figure  27. 


Figure  27.  Conceptual  sensitivity  field  using  1 800  UTC  23  January  2009 
analysis  temperature  as  the  analysis  state  variable  (. x, )  and  qvat 
0000  UTC  24  January  2009  as  the  forecast  metric  (J)  for  the  box 
centered  over  the  SLC  airport.  (After  Google  Maps) 

C.  SENSITIVITY  TESTING 

While  the  establishment  of  clear  sensitivity  relationships  is  important  for 
gaining  insight  into  the  relevant  atmospheric  dynamics,  this  information  alone 
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does  not  provide  sufficient  justification  for  placement  of  potential  atmospheric 
sensors.  In  order  for  an  additional  observation  to  be  of  use,  it  must  not  only 
exist  in  a  region  of  sensitivity,  it  must  also  be  shown  that  the  introduction  of  an 
observation  at  this  point  will  improve  a  forecast.  An  observation  added  at  a  point 
of  maximum  sensitivity  that  does  not  significantly  improve  the  forecast  is  of  little 
value. 

If  we  recall  the  equations  that  make  up  the  ESA  calculations,  maximum 
sensitivity  ( dJ/dx )  occurs  when  correlation  coefficient,  rxj,  is  strong  or  when 
spread  (standard  deviation)  of  analysis  state  variable,  crx,  is  small  for  a  fixed  crj. 
By  adding  an  observation  at  a  sensitivity  point,  we  expect  to  more  accurately 
sample  conditions  found  in  nature  (decreasing  ax),  thus  hopefully  improving  the 
analysis  that  follows  from  a  model.  Decreasing  crx  does  not  inherently  mean  a 
forecast  will  be  improved.  ESA  is  defined  by  both  rxj,  which  explains  the  linearity 
of  the  relationship  between  x  and  J,  and  the  ratio  of  crj  /  crx,  which  quantifies  the 
covariability  of  J  and  x.  The  contribution  to  sensitivity  by  these  variables  is  not 
equal.  Correlation  coefficient  can  only  serve  to  reinforce  sensitivity  values  when 
the  relation  between  xand  Jis  perfectly  linear  (r  =  +/-  1).  When  the  relationship 
is  not  perfectly  linear  (-1  <  r  <  +1),  the  correlation  coefficient  only  reduces  the 
overall  sensitivity  value.  To  maximize  sensitivity,  we  must  first  have  a  large  crj  / 
crx  ratio,  then  a  strong  associated  r-value.  Thus  in  choosing  to  introduce  an 
observation  at  a  position  of  max  sensitivity  defined  by  ESA,  we  cannot  be  certain 
that  an  improvement  in  forecast  accuracy  will  follow. 

1.  Perturbation  Tests 

Two  perturbation  tests  were  conducted  to  examine  whether  the 
introduction  of  an  observation  at  a  point  of  sensitivity  could  improve  a  later 
forecast  for  qv.  The  observation  is  represented  by  xs,  an  element  of  a  state 
vector  at  the  6-hour  ESA  most  sensitive  grid  point  X  =  121,  Y  =  52  (39. 5N, 
112.9W).  Analysis  ensemble  data  came  from  the  1800  UTC  23  January 
ensemble  mean  while  forecast  predictions  occur  at  0000  UTC  24  January. 
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a.  +1oT  Perturbation  Test 

The  first  test  introduced  a  perturbation  of  +1  a  temperature  for  xs  at 
the  first  /7-level.  This  test  simulated  the  introduction  of  a  real  observation.  The 
perturbation  was  introduced  by  increasing  the  ensemble  mean  temperature  value 
for  xs  by  1  standard  deviation  of  the  ensemble  mean  (+1oy  =  0.0516).  The  effect 
of  this  perturbation  was  spread  through  to  other  model  variables  in  the  horizontal 
and  vertical  through  simple  linear  regression  of  the  analysis  error  statistics 
sampled  by  the  analysis  ensemble.  This  linear  regression  ensured  a  balanced 
response  to  the  perturbation  and  subsequent  analysis. 

Figure  28  shows  the  1800  UTC  23  January  analysis  change  at  /7-1 
(in  K)  that  is  calculated  after  the  introduction  of  the  +1  oj  perturbation  at  the  most 
sensitive  grid  point  in  the  6-hour  ESA.  The  black  dot  indicates  the  location  where 
the  perturbation  was  introduced.  This  figure  shows  that  temperature  has 
increased  across  the  Salt  Lake  Valley,  and  at  KSLC.  While  it  may  be  expected 
to  see  the  strongest  response  to  the  perturbation  where  the  perturbation  was 
introduced  (xs),  because  no  covariance  localization  is  applied,  the  perturbation 
response  depends  entirely  on  the  covariance  of  the  background  temperature 
field.  It  has  already  been  discussed  how  significant  the  warm  air  advection  in 
and  around  KSLC  was  in  leading  to  the  dense  fog  at  KSLC  and  these  results 
further  support  that  reasoning.  Warmer  air  at  the  most  sensitive  grid  point  (+1  oy) 
helps  induce  a  stronger  inversion  over  the  southeast  GSL  as  the  warm  air 
overrides  the  cold  dome,  leading  to  a  later  dense  fog  intrusion  under  the 
strengthened  inversion.  Also  evident  in  Figure  28  is  a  region  of  strong  decrease 
in  temperatures  west  of  the  GSL.  This  area  represents  the  low  lying  Great  Salt 
Lake  Desert.  A  warm  perturbation  introduced  in  the  elevated  terrain  south  of  the 
Great  Salt  Lake  Desert  would,  with  southerly  flow,  likely  advect  warm  air  over  the 
PBL  in  this  region,  creating  or  strengthening  a  low  level  inversion  and  inducing 
colder  temperatures  than  in  surrounding  areas. 
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Figure  28.  Analysis  change  (in  K)  after  the  addition  of  +1  o  Q  perturbation 
applied  at  39.51  °N,  1 12.97 °W  at  the  first  /7-level.  The  perturbation 
point  is  indicated  by  the  black  dot.  Positive  values  indicate  an  increase  in  /7-1 
temperatures,  while  negative  values  indicate  a  decrease.  Note  the  strong 
positive  increase  in  Q  in  the  Salt  Lake  Valley,  directly  over  KSLC. 

Next  it  must  be  determined  what  change  in  forecast  is  created  by 

the  +1  oy  perturbation.  This  is  done  by  integrating  all  members  of  the  perturbed 

ensemble  model  forward  to  the  valid  forecast  time,  0000  UTC  24  January,  to  see 

if  these  new  forecasts  differs  from  the  unperturbed  forecasts.  The  focus  of  this 

test  is  to  examine  the  forecast  for  qv  at  KSLC  (J-box)  to  determine  what  change 

has  occurred.  Figure  29  shows  the  results  of  this  test.  The  graphic  represents 

the  difference  in  6-hour  forecasted  qv  (kg  kg"1)  between  the  forecast  created 

using  the  1800  UTC  23  January  +1oy  perturbed  analysis  and  the  forecast 

created  using  the  1800  UTC  23  January  unperturbed  analysis.  Note  that  this 

graphic  does  not  accurately  show  the  forecast  difference  within  the  entire  J-box 

as  it  only  displays  qv  at  the  lowest  model  level  and  the  J-box  value  includes  all 

points  over  the  first  two  model  levels. 
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Figure  29.  Forecast  difference  in  qv  (in  kg  kg"1)  at  0000  UTC  24  January  at  the 
lowest  model  level  (q. i)  following  the  introduction  of  a  +1  o  G 
perturbation  to  the  1800  UTC  23  January  ensemble  analysis  at  the 
most  sensitive  grid  point  (black  dot). 

The  total  difference  within  the  J-box  at  KSLC  between  the 
perturbed  and  unperturbed  forecast  is  difficult  to  quantify  by  visual  examination  of 
Figure  29.  The  J-box  appears  to  contain  positive  and  negative  areas,  while 
strong  positive  increases  (increased  qv)  are  evident  over  the  south  GSL  and 
south  of  the  J-box.  The  exact  forecast  difference  in  the  J-box  is  a  decrease  in 
(-7.3x1  O'6  kg  kg'1).  Though  it  would  be  preferred  to  see  an  increase  in  qv 
following  the  introduction  of  the  +1  oT perturbation,  there  is  a  clear  strong  positive 
response  across  the  southern  half  of  the  GSL  and  the  exact  change  at  KSLC  is 
only  weakly  negative.  The  known  forecast  timing  error  within  the  forecast  system 
(ending  the  fog  event  3-6  hours  early)  may  have  contributed  to  this  finding.  In 
addition,  the  significant  difference  in  and  around  the  J-box  suggests  that  a 
positive  forecast  change  exists  within  a  few  grid  points.  This  result  may  also 
indicate  the  presence  of  weak  nonlinearity  in  the  system. 
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Also  evident  in  Figure  29  are  several  areas  of  strong  increase  in  qv 
west  and  northwest  of  GSL,  along  the  border  of  northeast  Nevada  and  southern 
Idaho.  It  is  speculated  that  these  areas  are  associated  with  enhanced  vertical 
motion  along  significant  terrain  features  with  east-to-west  orientation  in  the 
region.  It  is  believed  that  the  warm  perturbation  induced  an  increase  in  southerly 
flow,  enhancing  convergence  along  the  elevated  terrain.  This  change  induced 
stronger  vertical  motion,  enhanced  adiabatic  cooling  and  condensation,  leading 
to  increases  in  qv concentrations. 

b.  “Truth”  Perturbation  Test 

A  different  procedure  was  then  completed  to  measure  the 
introduction  of  a  real  world  “truth”  observation  to  determine  if  a  more  accurate 
forecast  could  be  created.  Following  Evensen  (1994)  regarding  EnKF,  Anderson 
(2001)  regarding  Ensemble  Adjustment  Kalman  Filter  (EAKF),  and  Anderson  et 
al.  (2009)  regarding  DART  implementation  of  EnKF,  we  define  matrix  DA  for  our 
ensemble  as: 

=xf+K(y0-Hxf) 

where  xa  stands  for  the  analysis  background  variable,  T/the  forecast 

background  variable,  K  Kalman  gain,  yo  the  assimilated  observations,  and  H 
a  matrix  forward  operator. 

Matrix  DA  performed  with  a  single  observation  becomes 

xa  =  xf+K(y0-Hxf) 

with  yo  representing  the  single  observation  and  H  a  row  vector  forward 
operator.  K  is  now  defined  as 
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P  0  Pf  H 

HpoPfHT+(j2o 


where  p  represents  the  covariance  localization  matrix,  Pf  the  background  error 

variance  matrix,  °  the  Hadamard  or  Schur  product  between  the  two  arguments, 

2 

and  cr  0the  observation  error  covariance  matrix.  We  can  then  rewrite  our  DA  in 
scalar  form  with  single  scalar  additional  observation  as 


xa=xf  + 


P  cr 

2  Z  2  (?.-?/) 

p  cr  +  cr 

'  xv  XV  o 


2  2 

with  cr  ^representing  the  covariance  between  the  arguments,  cr  ^  the 

variance  between  the  arguments,  and  yf  the  forecast  value  of  yo,  (i.e.  HT Xf). 

This  equation  represents  how  the  additional  observation,  in  the  form  of  the 
perturbation,  is  assimilated,  allowing  the  contribution  of  one  single  observation  to 
be  measured. 

DA  for  this  perturbation  was  done  by  introducing  the  single 
observation  representing  the  truth  value  of  temperature  at  1800  UTC  23  January 
at  r/-i  for  xs.  The  single  observation  was  assimilated  into  the  96-member 
ensemble  as  defined  previously  in  this  section.  After  assimilating  the 
observation,  yn  =  xffor  truth,  the  model  background  crxof  9  at  xs  for  the  1800  UTC 
analysis  was  reduced  from  0.0516  to  0.0421.  The  effect  of  this  assimilated 
“truth”  observation  was  to  reduce  ox  and  thus  reduce  analysis  uncertainty. 

Figure  30  shows  the  difference  in  1800  UTC  23  January  analysis 
after  the  “truth”  observation  had  been  assimilated  at  the  most  sensitive  6-hour 
ESA  grid  point  (xs,  x  =  121,  y  =  52  (39. 5N,  112.9W)).  This  plot  is  obtained  by 
subtracting  the  analysis  of  the  unperturbed  model  run  from  the  analysis  of  the 
perturbed  model  run. 
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Figure  30.  1 800  UTC  23  January  analysis  change  (in  K)  after  the  addition  of 
“truth”  observation  applied  at  39.51  °N,  1 12.97°W  at  a].?.  The 
observation  point  is  indicated  by  the  black  dot.  Positive  values 
indicate  an  increase  in  A7.-1  temperatures,  while  negative  values 
indicate  a  decrease.  Note  the  very  strong  negative  response  in  6  in 
the  Salt  Lake  Valley,  directly  over  the  KSLC  J-box. 


Figure  30  shows  analysis  change  ranging  from  ±  20  K,  a  notable  increase 
in  magnitude  from  the  +1  oT  perturbation  which  yielded  analysis  change  ±  0.6  K. 
Importantly,  the  location  of  maximum  temperature  change  again  is  focused  on 
the  Salt  Lake  Valley,  in  and  around  KSLC,  though  the  analysis  change  in  this 
area  is  now  negative.  This  sign  change  is  expected,  as  the  assimilation  of  the 
truth  observation  {6x=  -0.0826  K)  lowered  the  background  temperature  at  xs(-1.1 
K).  The  similarity  in  geographic  position  of  analysis  change  between  +1oy 
perturbation  and  “truth”  observation  near  KSLC  is  important  to  our  results  and 
shows  that  the  ensemble  filter  is  a  modified  direct  linear  regression. 

Notably,  the  analysis  temperature  change  at  KSLC  after 
assimilating  the  “truth”  observation  was  significantly  greater  (-16.6  K)  than  the 
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analysis  change  at  xs  (-1.1  K).  The  significant  difference  in  analysis  change 
between  the  +1oy  perturbation  and  the  “truth”  perturbation  may  be  explained 
partially  by  the  difference  in  assimilation  procedures,  which  provided  different 
means  of  altering  initial  conditions.  The  more  significant  contribution  to  the 
difference  in  analysis  change  likely  comes  from  the  magnitude  of  the  difference 
between  the  +1  oT  perturbation  (+0.0516  K)  and  the  “truth”  observation  (-1.1  K). 
The  significantly  larger  “truth”  observation  leads  to  a  comparably  large  change  in 
analysis. 

Also  evident  is  an  area  of  strong  positive  analysis  change  west  of 
the  GSL.  This  region  represents  the  Great  Salt  Lake  Desert  that  was  identified 
with  the  +1  oT  perturbation  (Figure  28).  As  expected,  because  the  perturbation 
introduced  with  the  “truth”  test  was  negative,  this  decreased  the  amount  of  warm- 
air  advection  west  of  the  GSL,  thus  weakening  or  eroding  any  inversion  present 
in  this  area.  With  a  weakened  or  eroded  inversion,  cold  air  is  no  longer  trapped 
at  the  surface  as  efficiently  as  occurred  before  the  perturbation.  This  change 
results  in  an  increase  in  temperature  at  the  surface,  as  evidenced  by  the  positive 
analysis  temperature  change  in  this  area. 

As  predicted  by  our  ESA,  a  reduction  in  Q  at  xs  at  analysis  time 
(1800  UTC  23  January),  as  occurred  during  DA  of  the  “truth”  observation,  should 
predict  a  consequent  forecast  reduction  in  qvat  a  later  time  at  KSLC  (0000  UTC 
24  January).  Figure  31  shows  the  forecast  difference  that  occurs  between  the 
unperturbed  6-hour  forecast  and  the  6-hour  forecast  with  initial  conditions 
perturbed  by  the  “truth”  observation.  As  with  Figure  29,  this  graphic  represents 
only  the  lowest  model  level,  and  does  not  show  the  distribution  qv  difference  in 
the  entire  J-box. 


69 


-4 

x  10 


Figure  31 .  Forecast  difference  in  qv  (in  kg  kg"1)  at  0000  UTC  24  January  at  the 
lowest  model  level  {rj-i)  following  the  introduction  of  a  “truth” 
perturbation  to  the  1800  UTC  23  January  ensemble  analysis  at  the 
most  sensitive  grid  point  (black  dot) 

As  with  the  +1  oT perturbation  forecast  difference,  the  distribution  of 
the  forecast  qv  differential  across  the  J-box  appears  difficult  to  discern  by  visual 
inspection,  with  both  positive  and  negative  values  present.  Summation  of  qv 
within  the  J-box  reveals  an  overall  forecast  decrease  of  -4.765x1  O'4  kg  kg'1.  As 
expected,  the  magnitude  of  qv  forecast  difference  has  increased  by  an  order  of 
magnitude  from  the  +1  crT  perturbation,  roughly  matching  the  order  of  magnitude 
change  in  analysis  Q  at  xs.  In  addition,  the  concentration  of  higher  qv  that 
appeared  immediately  south  of  the  GSL  following  the  +1  oT  perturbation  (Figure 
29)  has  reversed  sign. 

Within  Figure  31,  a  concentrated  area  of  strong  (>  8.0x1 0"4  kg  kg'1) 
positive  forecast  change  is  evident  further  south  (-100  km)  of  the  GSL.  In  Figure 
16  {r\-i  ensemble  mean  qv),  this  same  region  appeared  as  a  pooling  area  for 
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water  vapor  between  0300  UTC  and  0600  UTC  24  January.  It  is  unclear  from 
Figure  31,  and  the  input  provided  by  the  previous  analyses  and  forecast 
difference  plots,  what  the  significance  of  this  region  is  and  how  it  relates  to  qvat 
KSLC.  It  is  speculated  that  a  decrease  in  warm-air  advection  in  this  region  may 
decrease  low-level  drying  that  would  have  resulted  from  the  advection  of  dry 
continental  air,  though  this  is  uncertain. 

These  results  are  positive  indications  that  the  introduction  of  an 
observation  produces  the  results  that  are  predicted  by  the  ESA  (e.g. 
increase/decrease  Q  at  xs  at  analysis  time  yields  increase/decrease  qv  at  later 
forecast  time).  The  next  step  is  to  analyze  how  these  results  compare  to  the 
theoretical  predictions  of  ESA. 

2.  ESA  Predicted  Sensitivity  Test 

ESA  calculations  provide  an  alternate  means  to  test  sensitivity  predictions. 
We  previously  defined  ESA  sensitivity  as  dJ/dx,  which  predicts  change  in  J  (6J) 
for  a  given  change  in  x  (6x).  Procedures  for  calculating  a  value  for  6x  differs 
based  upon  the  assimilation  methodology  used  which  is  defined  by  experiment 
purpose.  Once  the  procedure  for  6x  has  been  defined,  the  specific  value  of  6x 
may  be  multiplied  by  sensitivity  (6J  /6x),  the  product  of  which  is  an  expected 
change  in  J  (<5x),  our  sensitivity  prediction.  Figure  32  is  a  schematic  representing 
the  two  different  means  of  calculating  expected  <5Jfor  this  experiment. 
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Figure  32.  Schematic  representing  how  expected  change  in  J  (6J)  is  calculated 
for  (1)  +1cr  perturbation  and  (2)  “truth”  perturbation.  The  x-axis 
shows  temperature  (in  K)  at  the  most  sensitive  grid  point  in  the 
domain.  The  y-axis  shows  the  corresponding  forecast  qvat  KSLC  at 

a  later  time.  The  correlation  shown  is  idealized  as  perfectly  linear. 

a.  +1oT  Perturbation  Test 

To  simulate  an  observation,  8x  is  arbitrarily  defined  as  +J\o 
temperature  perturbation,  as  shown  in  Figure  32.  This  perturbation  moves  the 
ensemble  mean  9  at  the  most  sensitive  grid  point  by  +1cr.  Linear  regression  as 
previously  defined  was  performed  to  ensure  balanced  introduction  of  the 
perturbation.  With  6x  now  defined  as  +1  oT,  we  then  multiplied  this  value  (0.0564 
K)  by  the  6-hour  ESA  calculated  sensitivity  (6.61 53x1  O'4  kg  kg"1  K'1),  to  obtain  our 
expected  6J  as  3.3731  xl  O'5  kg  kg'1 . 

b.  “Truth”  Observation  Test 

As  previously  described,  our  matrix  DA  is  performed  as 

x“  =  xf  +K(y°  -Hxf) 

and  in  scalar  form  as 
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xa  -  xf  +  K(y°  —xf). 


This  can  be  then  rewritten  as: 

xa  —xf  =  K(y°  —xf) 

Q  J 

where  K(y  —x  )  represents  our  analysis  increment.  We  then  recognize 

^  a  ^  a 

/~j~ v*  /]' v- 

xa  -xJ  as  — .  Multiplying  —times  our  previously  defined  sensitivity  value 

dy  dy 

function,-^-,  we  return 

dxa 

8xa  8J  _8J 
dy  8xa  8y 

dJ 

which  provides,  — ,  an  expected  change  in  forecast  variable  (J),  given  a  change 

dy 

in  observation  (y).  The  input  iorxa  -xf ,  or  6x,  is  provided  as  the  ©-difference  at 
rj-i  between  the  1800  UTC  23  January  “truth”  temperature  (Figure  32  -  Ttruth),  and 
the  1800  UTC  23  January  analysis  mean  background  field  (Figure  32  -  Tens)  at 

dxa 

the  most  sensitive  grid  point.  This  value  (  —  =  -1.0528  K)  is  then  multiplied  by 

dy 

dJ 

the  sensitivity  value  at  the  6-hour  ESA  most  sensitive  grid  point  ( —  = 

dx“ 

d  I 

6.6153x10'4  kg  kg"1  K'1).  The  resulting  value  ( —  =  -6.9646x10"4  kg  kg'1)  is  the 

dy 

expected  change  in  6J  (qv)  at  KSLC  at  0000  UTC  24  January. 

3.  Actual  Versus  Expected  Forecast  Change 

If  we  consider  the  simulations  performed  through  perturbation 
testing  (section  III.C.1.  and  III.C.2.)  as  the  actual  change  in  forecast  and  then 
consider  the  ESA  predicted  test  as  the  expected  change  in  forecast,  we  can 
compare  these  results  to  determine  the  validity  of  the  ESA. 
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+1  a  Perturbation 
(K) 

"Truth" 

Observation  (K) 

Perturbation  Test 

(Actual  Change) 

-7.30x1  O'6  kg  kg'1 

-4.77x1  O'4  kg  kg'1 

ESA  Prediction  Test 

(Expected  Change) 

3.73x1  O'5  kg  kg'1 

-6.97x1  O'4  kg  kg'1 

Table  3.  Comparison  of  Actual  versus  Expected  qv  changes  in  kg  kg"1  at 
KSLC  at  0000  UTC  24  January  with  (1 )  +1  o  Perturbation  (K)  and 
with  (2)  “Truth”  Observation  (K)  at  6-hour  ESA  most  sensitive  grid 
point  at  1800  UTC  23  January. 

Table  3  shows  the  results  of  each  test.  First  looking  at  the  +1cr 
perturbation,  it  is  evident  that  the  actual  change  (-7.30x1  O'6  kg  kg'1)  and  the 
expected  change  (3.73x1  O'5  kg  kg'1)  differ  in  both  sign  and  magnitude.  It  is 
speculated  that  poorly  correlated  values  presented  by  the  actual  perturbation  test 
can  be  attributed  to  deficiencies  in  how  the  observation  is  assimilated.  As 
described  earlier,  the  perturbation  is  introduced  by  increasing  the  ensemble 
mean  9  value  +1  o,  followed  by  linear  regression  to  distribute  the  new  value.  This 
process  may  not  fully  capture  the  impact  of  the  perturbation,  thus  introducing 
errors  in  the  system  which  contributes  to  error  in  forecast  values. 

Assimilating  the  “truth”  observation  using  DA  yields  amore 
encouraging  results  with  agreement  in  sign  (negative)  and  magnitude  (xIO  4  kg 
kg"1).  Because  our  “truth”  observation  had  reduced  our  ensemble  mean  6  for  xs, 
a  subsequent  decrease  in  qvat  KSLC  at  a  later  time  is  expected,  and  is  shown  to 
be  true  in  Table  3.  This  agreement  is  significant  and  informs  us  that  our  ESA  can 
successfully  predict  changes  in  fog  at  KSLC. 

D.  ENSEMBLE  SIZE  TESTING 

The  establishment  of  the  size  of  the  ensemble  system  (member  of 
ensemble  members)  needed  to  perform  a  reliable  ESA  is  one  of  the  primary 
focuses  of  this  study.  Identifying  this  information  helps  to  establish  the 
computational  requirements  necessary  to  implement  a  reliable  ESA.  This  factor 


74 


is  important  when  determining  whether  or  not  an  NWP  computing  center  can 
operationally  implement  such  measures  into  their  DA  process.  Even  if  not  used 
operationally,  the  need  to  use  the  appropriate  number  of  ensemble  members  is 
critical  to  building  an  accurate  sensitivity  study  for  placement  of  future 
observations. 

Chilcoat  (2012)  identified  the  requirement  for  an  ensemble  system  to 
consist  of  30-40  members  to  skillfully  determine  sensitivity  of  forecast  fog  at 
KSLC.  In  addition  to  the  use  of  a  different  initial  condition  parameter,  Chilcoat’s 
research  did  not  implement  the  stringent  requirement  for  Gaussianity  that  was 
implemented  for  this  research;  therefore  different  results  are  to  be  expected. 

Ensemble  sub-selection  procedures  were  detailed  in  section  III.  Results 
were  obtained  using  ensemble  sizes  of  96,  80,  64,  and  48  members  for  both  6- 
hour  and  12-hour  sensitivity  analyses.  Sensitivities  {dJ/dxj)  were  calculated  using 
initial  condition  T  (K)  as  x,  and  forecast  qv  (kg  kg"1)  for  Jthe  12km  x  12  km  J-box 
surrounding  KSLC.  Comparison  of  sensitivity  results  between  ensemble  systems 
occurred  through  a  visual  analysis  of  sensitivity  plots  and  numerical  analysis  of 
sensitivity  values.  Continuity  between  ensemble  systems  was  established  by 
identifying  similar  geographical  positioning  and  magnitude  of  +/-  sensitivities 
across  the  domain.  The  10  strongest  sensitivities  were  identified,  by  grid  point, 
for  each  size  ensemble  system  at  both  6-hour  and  12-hour  times.  The  three 
strongest  sensitivity  points  are  identified  in  Figures  33  and  34.  All  sensitivities 
failing  the  95%  confidence  interval  test  are  masked. 

1.  6-hour  ESA 

Figure  33  shows  the  result  of  the  ESA  for  all  four  6-hour  ensemble 
systems.  The  96-member  ensemble  system  served  as  the  standard  which  the 
other  systems  were  compared  against.  Features  of  the  96-member  system  were 
defined  in  section  IV. B.  The  most  sensitive  point  in  the  96-member  ensemble 
exists  at  (39.5072N,  112.9197W)  and  is  identified  at  Point  1,  96-member  plot, 
Figure  18.  Notable  in  the  96-member  ensemble  is  strong  agreement  within  the 
model  system  of  the  location  of  greatest  sensitivity,  with  seven  of  the  top  ten 

sensitivities  existing  in  and  around  the  most  sensitive  grid  point. 
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Figure  33.  Plot  of  6-hour  sensitivity  ( dJ/dx, )  for  96,  80,  64,  and  48-member 
ensemble  systems  using  initial  condition  T (K)  from  1800  UTC  23 
January  2009  and  forecast  qv  (kg  kg"1)  at  KSLC  at  0000  UTC  24 
January  2009.  Units  are  kg  kg'1  K"1.  Number  boxes  identify  grid 
points  of  the  strongest  three  sensitivities  in  each  plot.  Markers 
indicating  (+)  or  (-)  represents  a  positive  or  negative  sensitivity  that 
exists  outside  the  graphic.  Top  left  number  represents  number  of 
ensemble  members  used  to  perform  ESA.  Note:  scales  are  not 

normalized. 


The  80-member  system  agreed  relatively  well  with  positioning  of  96-hour 
features.  Areas  of  +/-  sensitivity  retain  similar  geographical  arrangement  which 
suggest  good  agreement  between  the  ensembles.  Magnitude  of  the  features 
does  not  match  as  well,  with  only  the  strongest  sensitivity  grid  point  retaining  its 
position,  and  overall  reduction  in  sensitivity  values  with  the  most  sensitive  grid 
point  reducing  in  strength  from  6.61 5x1  O'4  kg  kg"1  K"1  to  3.021x1  O'4  kg  kg'1  K'1. 
Note  also  the  reduction  in  scale  between  96-  and  80-member  ensemble  systems. 
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In  addition,  the  next  two  strongest  grid  points  existed  in  different  areas  of  the 
map,  indicating  a  breakdown  in  the  agreement  of  the  most  sensitive  region  that 
existed  in  the  96-member  system. 

The  64-member  ensemble  also  retained  a  similar  geographic  arrangement 
as  the  96-hour  member  with  the  most  salient  features  still  present.  This  system 
begins  to  experience  more  magnitude  displacement  as  the  three  most  sensitive 
grid  points  are  well  southeast  of  the  GSL,  though  they  occur  in  relative  proximity 
to  each  other.  The  location  of  maximum  sensitivity  in  the  96-member  system 
now  ranks  ninth  in  magnitude  with  sensitivity  of  4.11 4x1 0"4  kg  kg"1  K'1.  In 
addition  some  different  +/-  groupings  begin  to  emerge  that  differ  from  those  in  the 
96-  and  80-member  systems.  Because  of  the  deterioration  in  magnitude  and 
position  of  sensitivities,  this  is  likely  below  the  minimum  size  ensemble  needed  to 
successfully  implement  a  6-hour  ESA. 

The  48-member  ensemble  also  retains  the  most  salient  sensitivities  from 
the  96-hour  member,  though  most  of  the  sensitivities  fail  the  Cl  test. 
Interestingly,  the  response  at  the  96-member  region  of  maximum  sensitivity 
improves  from  the  64-member  ensemble,  as  the  second  and  ninth  most  sensitive 
grid  points  emerge  in  this  area.  The  system  has  lost  little  sensitivity  magnitude, 
as  the  sensitivity  at  Point  2  (Figure  33  -  bottom  right)  is  4.239x1 0"4  kg  kg'1  K"1, 
which  is  second  in  magnitude  to  the  96-member  ensemble  responses  at  this 
point.  This  seemingly  counterintuitive  increase  in  sensitivity  can  be  attributed  to 
sampling  error.  As  ensemble  size  decreases  covariances  can  become 
overestimated,  while  variances  can  become  underestimated.  This  occurs  as  the 
increasingly  smaller  ensemble  samples  a  smaller  portion  of  the  true  population, 
drawing  a  less  accurate  representation  of  the  system  parameters  and  increasing 
sampling  error. 

2.  12-hour  ESA 

The  overall  response  of  the  12-hour  ensemble  sub-systems  (Figure  34) 
varied  significantly  from  the  6-hour.  The  96-member  ensemble  again  served  as 
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the  control  “best-response”  for  comparison.  Explanation  of  the  12-hour  96- 
member  ensemble  sensitivity  is  found  in  section  IV.B.  The  primary  features  of 
the  96-member  12-hour  system  indicate  strong  positive  sensitivities  in  the 
Central  Utah  Desert  Range  and  in  lower  lying  valleys  south  and  east  of  the  GSL. 
The  point  of  maximum  sensitivity  exists  in  the  same  general  area  as  the  6-hour 
system  (identified  as  points  1,  2,  and  3,  Figure  34  -  top  left),  with  maximum 
sensitivity  of  4.1 92x1 0"4  kg  kg"1  K"1. 


Figure  34.  Plot  of  1 2-hour  sensitivity  ( dJ/dx, )  for  96,  80,  64,  and  48-member 
ensemble  systems  using  initial  condition  T (K)  from  1800  UTC  23 
January  2009  and  forecast  qv  (kg  kg'1)  at  KSLC  at  0000  UTC  24 
January  2009.  Units  are  kg  kg"1  K'1.  Number  boxes  identify  grid 
points  of  the  strongest  three  sensitivities  in  each  plot.  Markers 
indicating  (+)  or  (-)  represent  a  positive  or  negative  sensitivity  that 
exists  outside  the  graphic.  Top  left  number  represents  number  of 
ensemble  members  used  to  perform  ESA.  Note:  scales  are  not 

normalized. 
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The  80-member  ensemble  system  exhibited  notable  difference  from  the 
96-member  system.  While  the  region  of  positive  sensitivities  in  the  Central  Utah 
Desert  Ranges  was  evident,  this  region  produced  only  the  sixth  strongest 
sensitivity  in  the  domain.  The  strongest  sensitivities  emerged  on  the  Front 
Range  of  the  Wasatch  Mountains,  immediately  east  of  KSLC.  Values  here 
exhibited  maximum  sensitivity  of  3.964x1  O'4  kg  kg"1  K"1.  This  system  also 
produced  a  number  of  seemingly  spurious  +/-  sensitivity  regions,  when  compared 
against  the  96-member  system.  This  result  is  to  be  expected,  when  considering 
the  mesoscale  conditions  that  existed  over  the  12-hour  period.  While  warm  air 
advection  was  predominant  across  central  and  southern  Utah,  northern  Utah 
experienced  greater  variability  in  forcing  as  the  effect  of  the  mid-level  disturbance 
was  replaced  by  southerly  flow. 

The  sensitivity  profile  of  the  64-member  ensemble  did  not  resemble  the 
96-member  system  in  geographical  area  or  magnitude.  No  previously 
recognized  sensitivity  features  were  evident  with  the  64-member  system.  Of 
significance  was  the  magnitude  that  did  exist  with  the  system.  The  most 
sensitive  point,  well  southeast  of  GSL  and  identified  as  Point  1 ,  Figure  34 
(bottom  left),  had  sensitivity  of  7.276x1 0"4  kg  kg'1  K'1,  which  was  the  strongest 
value  exhibited  at  any  point  during  this  study.  The  lack  of  coherent  sensitivity  in 
this  system  indicates  a  significant  discontinuity  with  the  96-  and  80-member 
systems. 

The  48-member  ensemble  system  exhibited  even  less  coherence  in 
sensitivity  response,  with  most  points  on  the  domain  failing  the  Cl  test.  Though 
the  strongest  sensitivity  value  is  comparable  to  other  systems  at  4.344x1  O'4  kg 
kg"1  K'1,  this  system,  like  the  64-member  system,  is  not  likely  useful  as  an 
accurate  ESA  system.  Further  examination  will  explore  this  proposal. 

3.  Linearity  of  the  6-hour  Positive  Sensitivities 

Figure  35  is  the  scatterplots  for  the  6-hour  ESA  created  with  96-,  80-,  64-, 
and  48-member  ensemble  systems.  The  x-axis  in  each  scatterplot  in  Figure  35 
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represents  the  ensemble  mean  analysis  potential  temperature  (in  K)  at  1800 
UTC  23  January  2009  at  the  lowest  model  level  at  the  most  sensitive  grid  point 
(39.5072N,  112.9197W)  in  the  6-hour  ESA.  The  y-axis  in  each  scatterplot 
represents  the  forecasted  qv  (in  kg  kg-1)  in  the  J-box  surrounding  KSLC  at  0000 
UTC  23  January  2009. 
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Figure  35.  Scatterplots  of  analysis  q.-i  temperature  and  forecasted  qv  used  in 
the  ESA  as  ensemble  size  is  reduced  from  96  to  48  members.  The  x- 
axis  represents  the  1800  UTC  23  January  analysis  temperature  (in 
K)  at  the  most  sensitive  grid  point  (39.5 °N,  1 1 2.9 °W).  The  y-axis 
represents  the  forecasted  qv  value  (in  kg  kg"1)  six  hours  later,  valid 
0000  UTC  24  January.  Also  plotted  are  correlation  coefficient  (r)  and 

sensitivity  ( dJ/dx ). 


The  top  left  panel  in  Figure  35  represents  the  96-member  6-hour 
scatterplot  which  was  discussed  in  section  IV.3.B.  The  top  right  panel  in  Figure 
35  is  the  scatterplot  for  the  80-member  ensemble.  This  panel  indicates  that  the 
80-member  ensemble  has  weaker  sensitivity  than  the  96-member  ensemble  (a 
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reduction  from  6.61 53x1 0"4  kg  kg"1  K'1  to  1.7780x1  O'4  kg  kg'1  K"1).  The  80 
member  ensemble  also  exhibits  less  linearity  than  the  96-member  ensemble  (a 
reduction  in  correlation  from  r  =  0.55  to  r  =  0.22).  The  bottom  two  panels  in 
Figure  35  represent  the  64  and  48-member  ensemble,  respectively.  Each  of 
these  scatterplots  shows  a  significant  reduction  in  sensitivity  (an  order  of 
magnitude  reduction)  and  very  poor  linear  correlation  between  Q  and  qv  (r-values 
<  0.1).  This  significant  reduction  in  sensitivity  and  linearity  between  80  and  64 
ensemble  members  appears  to  define  the  lower  limit  of  ensemble  size  that  is 
needed  to  perform  a  6-hour  ESA  for  this  case. 

Figure  36  provides  a  depiction  of  both  correlation  and  sensitivity  for  the  6- 
hour  ESA  for  each  ensemble  member.  In  addition,  the  green  bars  in  Figure  36 
indicate  the  upper  and  lower  boundaries  of  the  95%  confidence  interval  for  the 
correlation  coefficient  of  each  member.  These  boundaries  represent  the 
sampling  error  associated  with  each  ensemble  member.  Reduction  in  ensemble 


size  is  known  to  correspond  with  increasing  sampling  error. 


Figure  36.  6-hour  ESA  chart  of  correlation  in  red  and  sensitivity  in  blue  versus 
the  number  of  ensemble  members  for  a  single  point  of  positive 
sensitivities  at  39.5 °N,  1 1 2.9 °W.  Green  bars  indicate  the  high  and 
low  boundaries  of  the  95%  confidence  interval  and  represent  the 
associated  sampling  error. 
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Figure  36  confirms  the  results  shown  in  Figure  35.  As  expected, 
ensembles  smaller  than  80  members  exhibit  significant  reduction  in  sensitivity 
and  linearity.  In  addition,  sampling  error  associated  with  the  r-value  of  each 
ensemble  member  increases,  also  as  expected.  Of  importance  is  the  error 
associated  with  the  80-member  ensemble.  The  green  error  bars  indicate 
significant  variability  (r  =  +/-  0.20)  in  the  correlation  coefficients  associated  with 
this  ensemble.  This  result  means  there  is  high  variability  in  the  linearity  that 
accompanies  the  80  member  ensemble  and  that  caution  should  be  taken  when 
implementing  an  ESA  with  fewer  than  96-members.  To  achieve  the  most  robust 
and  reliable  results  for  the  6-hour  ESA,  a  96-member  ensemble  is  preferred. 

4.  Linearity  of  the  12-hour  Positive  Sensitivities 

Previous  results  presented  in  section  IV.3.B.  showed  the  significant 
reduction  in  linearity  that  occurred  when  assessing  the  sensitivity  of  the  12-hour 
ESA  at  the  most  sensitive  6-hour  ESA  grid  point.  Figure  37  contains  the 
scatterplots  that  accompany  ensemble  sub-selection  of  this  grid  point.  Linear 
associations  range  from  weak  and  positive  (r  =  0.10,  64  members)  to  lacking 
linearity  (r  =  0.023,  96  members;  r  =  0.025,  48  members),  to  moderate  and 
negative  (r  =  -0.3176,  80  members).  This  data  further  emphasizes  the  significant 
variability  that  accompanies  the  sensitivity  and  the  linearity  between  6  and  12- 
hour  time  frames.  Though  this  data  point  was  masked  for  the  12-hour  ESA  for 
failing  the  95%  Cl  test  (P-value  =  0.8218),  this  variability  indicates  the  limitations 
of  sensitivity  data  from  a  single  grid  point,  especially  when  ESA  is  performed  with 
less  than  96  members. 
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Figure  37.  Scatterplots  of  analysis  q.-\  temperature  and  forecasted  qv  used  in 
the  ESA  as  ensemble  size  is  reduced  from  96  to  48  members.  The  x- 
axis  represents  the  1800  UTC  23  January  analysis  temperature  (in 
K)  at  the  12-hour  ESA  indicated  most  sensitive  grid  point  (39.5°N, 

1 1 2.9 °W).  The  y-axis  represents  the  forecasted  qv  value  (in  kg  kg"1) 
12  hours  later,  valid  0000  UTC  24  January.  Also  plotted  are 
correlation  coefficient  (r)  and  sensitivity  ( dJ/dx ). 


Also  discussed  in  IV.3.B.  was  the  successful  increase  in  the  12-hour 
linearity  that  occurred  following  a  change  from  the  6-hour  ESA  most  sensitive 
grid  point  to  the  12-hour  ESA  most  sensitive,  and  second  most  sensitive  grid 
points.  Figures  38-41  show  the  accompanying  changes  in  sensitivity  and 
linearity  that  occur  during  ensemble  sub-selection.  Figure  38  provides 
scatterplots  for  all  ensemble  four  ensemble  systems  comparing  1200  UTC  23 
January  analysis  potential  temperature  at  the  12-hour  ESA  most  sensitive  grid 
point  and  0000  UTC  24  January  qvat  KSLC.  Figure  39  presents  the  r-values 
(with  associated  sampling  error)  and  sensitivities  for  this  same  data  set.  Figure 
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40  and  Figure  41  provide  equivalent  graphs  for  the  12-hour  ESA  second  most 
sensitive  grid  point. 

These  figures  indicate  that  sub-selection  of  ensembles  results  in  less 
linearity  with  decreasing  ensemble  size,  and  less  overall  linearity  when  compared 
with  equivalent  6-hour  ESA  sub-systems.  R-values  for  80-member  ensembles 
fluctuate  from  lacking  linear  relationships  (r  =  0.087,  most  sensitive  12-hour  ESA 
grid  point)  to  weakly  linear  (r  =  0.253,  second  most  sensitive  12-hour  ESA  grid 
point).  64  and  48-member  ensembles  for  both  top  12-hour  sensitivity  points  lack 


Figure  38.  Scatterplots  of  analysis  /7-1  temperature  and  forecasted  qv  used  in 
the  ESA  as  ensemble  size  is  reduced  from  96  to  48  members.  The  x- 
axis  represents  the  1200  UTC  23  January  analysis  temperature  (in 
K)  at  the  12-hour  ESA  indicated  most  sensitive  grid  point  (8km  N  and 
4  km  E  of  39.5 °N,  1 1 2.9 °W).  The  y-axis  represents  the  forecasted  qv 
value  (in  kg  kg"1)  twelve  hours  later,  valid  0000  UTC  24  January. 
Also  plotted  are  correlation  coefficient  (r)  and  sensitivity  ( dJ/dx ). 
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12-hr  Sensitivities  8km  N,  4km  E  of  39.5°N,  112.9°W 
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Figure  39.  12-hour  ESA  chart  of  correlation  in  red  and  sensitivity  in  blue  versus 
the  number  of  ensemble  members  for  a  single  point  of  positive 
sensitivities  at  the  most  sensitive  grid  point  for  the  1 2-hour  ESA,  8km 
N  and  4km  E  of  39.5 °N,  11 2.9 °W.  Green  bars  indicate  the  high  and 
low  boundaries  of  the  95%  confidence  interval  and  represent  the 
associated  sampling  error. 
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Figure  40.  Same  as  Figure  38,  with  12-hour  ESA  indicated  second  most 
sensitive  grid  point  (28km  N  and  1 2  km  E  of  39.5 °N,  112.9 °W) 
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Figure  41.  Same  as  Figure  39,  for  a  grid  point  28  km  N  and  12  km  E  of  39.5°N, 

1 1 2.9 °W.  This  grid  point  is  the  second  most  sensitive  grid  point  in 

the  12-hour  ESA. 

5.  Ensemble  Size  Conclusions 

The  results  of  the  6-hour  and  12-hour  ESA  ensemble  sub-selection  testing 
indicate  the  need  for  a  96-member  ensemble  in  order  to  produce  the  most 
statistically  significant  and  linear  sensitivities.  The  64-member  and  48-member 
ensembles  provided  sensitivity  results  too  weak  to  utilize  or  with  insufficient 
linearity  to  prove  useful  for  ESA.  The  80-member  ensembles  exhibited  weak 
linear  correlations,  but  also  with  weaker  overall  sensitivity  response  and  with 
greater  error  associated  with  the  correlation.  The  variability  exhibited  between 
grid  points  for  80-member  ensembles  indicates  that  correlation  derived  from  an 
ensemble  of  this  size  is  subject  to  significant  error  variation  (12-hour  r-value  error 
=  ±0.21-0.22).  Though  not  significantly  greater  in  magnitude  than  the  96- 
member  ensemble  (12-hour  r-value  error  =  ±0.17-0.19),  this  error  becomes  more 
critical  with  80-ensemble  members  due  to  the  low  correlation  values  returned. 
This  variation  at  low  correlation  values  means  that  the  ESA  calculation  returned 
with  an  80-member  ensemble  (12-hour  ESA  r-value  =  0.09-0.25)  may  be 
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moderately  linear  or  may  not  contain  linearity.  96-member  ensemble  sensitivities 
by  contrast,  which  show  greater  linearity  (12-hour  ESA  r-value  =  0.31-0.39)  and 
less  correlation  error,  provide  a  more  accurate  and  reliable  sensitivity  response. 
If  an  80-member  ensemble  is  to  be  used  to  perform  ESA,  care  must  be  taken  to 
ensure  that  non-linearities  are  not  an  issue. 

The  significant  variation  in  sensitivity  and  linearity  over  6-hour  and  12-hour 
time  frames  for  the  all  ensemble  systems  considered  provides  further  support  for 
the  need  for  a  spatially  expansive  network  when  considering  the  placement  of 
additional  observational  sensors.  The  single  source  of  failure  potential  that 
accompanies  the  data  from  a  single  grid  point  is  evident  with  these  results  and 
underscores  the  complexity  of  the  statistical  and  dynamic  relationships  in  this 
mesoscale  complex  terrain. 
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V.  CONCLUSIONS  AND  FUTURE  WORK 


This  study  provided  a  second  investigation  of  the  efficacy  of  ESA  in 
complex  terrain,  at  high  model  resolution,  and  in  the  mesoscale.  The  terrain 
involved  the  Great  Salt  Lake  Basin,  a  region  of  dynamic  changes  in  elevation  and 
hydrologic  composition  in  northern  Utah.  The  high  resolution  model  used  in  this 
study  was  the  WRF-ARW,  implemented  at  4-km  grid  spacing  on  the  innermost 
domain  of  a  nested  grid.  A  96-member  ensemble  of  this  model  was  then 
created.  ESA  was  performed  for  a  single  dense  fog  event  that  occurred  at  Salt 
Lake  City  International  Airport  (KSLC)  in  the  overnight  hours  of  23-24  January 
2009.  Our  analysis  of  the  ESA  indicated  a  strong  and  dynamically  relevant 
connection  between  observed  low  level  temperatures  and  forecasted  fog  at  the 
airport. 

Sensitivity  testing  began  by  performing  an  ESA  using  multiple  different 
initial  condition  parameters  at  1800  UTC  23  January  and  forecast  parameters  at 
0000  UTC  24  January.  Initial  results  indicated  that  temperature,  specifically 
potential  temperature  (0)  at  the  lowest  model  level  {rj.i),  provided  the  clearest 
and  strongest  6-hour  sensitivities  to  forecast  fog  (qv)  at  KSLC.  Other  forecast 
parameters,  including  qv,  total  dry  air  mass  in  a  column,  and  u-  and  v-  wind 
components,  failed  to  provide  results  as  robust  as  0.  The  ESA  temporal  scale 
was  then  increased  to  12  hours  to  determine  the  applicability  of  these  results  for 
longer  time  periods.  Using  1200  UTC  23  January  0  and  0000  UTC  24  January 
qv,  it  was  evident  in  the  12-hour  ESA  that  0  continued  to  provide  a  generally 
spatially  consistent,  yet  weaker  sensitivity  response  than  during  the  initial  6-hour 
test.  Following  these  tests,  it  was  concluded  that  0  would  be  the  most  reliable 
predictor  variable  to  use  in  performing  this  ESA. 

To  meet  the  Gaussian  distribution  requirements  inherent  in  ESA  as  a 
system  based  on  linear  regression,  we  then  tested  our  data  sets  to  ensure  they 
were  normally  distributed.  An  element  of  a  state  vector  at  a  shared  single  most 
sensitive  grid  point  (xs)  was  identified  from  6-hour  and  12-hour  ESA  results,  and 
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distribution  of  rj-i  initial  condition  Q  (xi)  for  each  96-member  ensemble  was 
analyzed  for  Gaussianity.  A  Lilliefors  test  was  implemented  to  perform 
Gaussianity  testing  and  results  indicated  that  the  6-hour  ESA  Xi  ensemble 
distribution  was  Gaussian.  The  12-hour  Lilliefors  test  indicated  that  the  x, 
ensemble  distribution  at  this  same  point  was  non-Gaussian,  though  several 
nearby  grid  points  showed  valid  Gaussianity.  It  was  concluded  that  positive  6- 
hour  Gaussianity  results,  combined  with  positive  12-hour  Gaussianity  at  close- 
proximity  grid  points,  was  sufficient  to  justify  our  input  parameter  distribution  as 
normal.  The  same  test  was  then  performed  for  the  forecast  parameter  (J)  at 
KSLC.  Rather  than  measuring  Gaussianity  at  a  single  grid  point  as  was  done  for 
it  was  necessary  to  measure  Jin  a  12-km  x  12-km  box  (4x4  grid  spacing)  over 
the  lowest  two  model  levels,  in  order  to  find  a  normal  distribution  of  J.  This  area- 
averaged  J  provided  more  robust  results  by  averaging  out  sampling  error  and 
noise. 


Once  Gaussianity  testing  confirmed  our  forecast  distributions  were 
Gaussian,  we  moved  forward  with  conducting  the  6  and  12-hour  ESA  tests  under 
these  new  conditions.  In  order  to  account  for  statistical  sampling  error,  this  data 
was  subjected  to  a  confidence  interval  test  at  the  95%  level.  Sensitivities  failing 
the  confidence  interval  test  (p>0.05)  were  removed  from  the  analysis.  The 
results  of  the  6-hour  ESA  indicated  a  response  closely  related  to  warm  air 
advection.  Strong  positive  and  negative  sensitivities  to  temperature  were  found 
at  geographic  areas  likely  to  see  significant  influence  of  warm  air  advection.  It 
was  believed  this  suggested  a  relation  to  strengthening  of  a  low  level  inversion 
which  induced  fog  at  KSLC  at  a  later  time.  Results  of  the  12-hour  ESA  showed 
an  overall  weaker  sensitivity  response,  but  one  still  influenced  by  warm  air 
advection.  These  results  showed  additional  influence  from  the  mid-level 
disturbance  that  passed  through  northern  Utah  on  23  January.  This  disturbance 
noticeably  alters  the  distribution  of  sensitivities,  though  the  dynamic  signal  of 
both  forcing  mechanisms,  warm  air  advection  and  mid-level  disturbance,  can  be 
distinguished.  The  non-linearities  in  u-  and  v-  component  wind  introduced  by  the 
mid-level  disturbance  had  prompted  Chilcoat  (2012)  to  move  the  focus  of  his 
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study  6  hours  later,  to  avoid  complication.  This  study  did  not  make  such 
changes,  as  9  provided  a  sufficiently  robust  result. 

Examination  of  the  linearity  of  the  data  sets  revealed  that  6-hour 
correlation  between  x,  and  J  existed  at  a  moderate  level  (r  =  0.55).  The  12-hour 
ESA  indicated  that  the  linearity  broke  down  when  the  6-hour  ESA  xs  was  used, 
but  could  be  restored  to  weak  or  moderate  linearity  (r  =  0.33  -  0.39)  by  moving  xs 
to  12-hour  ESA  most  sensitive  grid  points.  This  finding  shows  that  under 
influence  of  complex  mesoscale  forcing  over  varying  temporal  scales,  sensitivity 
locations  exhibit  dynamic  spatial  variation.  This  implies  that  a  future  observing 
network  must  include  a  mobile  observing  platform,  or  be  distributed  with  sufficient 
density  at  known  sensitivity  regions  as  to  adequately  observe  changes  in  initial 
conditions  known  to  relate  to  changes  in  forecast  model  output. 

Tests  were  then  conducted  to  determine  whether  the  addition  of  an 
observation  y0  could  accurately  predict  a  change  in  J  at  a  later  time.  Tests  were 
performed  with  a  +1cr  perturbation,  to  simulate  an  observation,  and  with  a  “truth” 
perturbation  to  represent  the  effects  of  assimilating  an  observation.  ESA  also 
allows  for  sensitivity  testing  using  ensemble  statistics  and  “truth”  data  to  create 
“expected”  forecast  changes,  done  without  actually  introducing  a  perturbation. 
Comparison  of  “truth”  sensitivity  testing  indicates  the  introduction  of  a  6 
observation  at  rj-i  for  xs  will  accurately  predict  change  in  qv  6-hours  later  at  KSLC. 
This  result  is  significant  in  that  it  shows  potential  for  ESA  to  be  used  to  place  an 
observation (s)  which  can  then  improve  forecasts  for  parameters  of  specific 
interests. 

The  final  assessment  performed  was  to  analyze  how  ensemble  sub¬ 
sampling  affected  sampling  error,  the  results  of  which  could  identify  the  minimum 
size  ensemble  needed  to  perform  ESA.  Both  6  and  12-hour  ESA  subsampling 
was  conducted  and  both  indicated  a  strong  preference  for  an  ensemble  of  96- 
members  for  best  results.  While  80-member  ensembles  could  potentially  be 
used  to  achieve  similar  results,  the  weakness  of  the  sensitivity  and  the  amount  of 
increased  error  associated  with  the  correlations  are  significant  shortcomings  for 
an  ensemble  of  this  size. 
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A.  COMPARISON  TO  CHILCOAT  (201 2)  RESULTS 


Qualitative  comparison  between  this  study  and  Chilcoat  (2012)  indicates  a 
general  agreement  in  the  evolution  of  the  dense  fog  event  at  KSLC.  Both  studies 
showed  similar  forecasting  accuracy  for  the  density  of  qv  at  KSLC  and  both 
indicated  a  tendency  to  begin  and  end  the  dense  fog  event  3-6  hours 
prematurely.  Differences  in  conclusions  between  the  two  studies  primarily  arise 
because  of  the  difference  in  ESA  forecast  time  frames.  Because  of  the  time  shift 
to  a  later  hour,  Chilcoat’s  u-  and  v-wind  ESA  results  and  resulting  conceptual 
model  identify  the  locations  and  dynamic  features  most  responsible  for  advecting 
the  fog  into  KSLC.  This  study,  in  contrast,  identified  the  initial  conditions  most 
responsible  for  the  atmospheric  profile  that  made  the  creation  of  the  dense  fog 
event  possible.  Together,  these  studies  provide  a  thorough  conceptual 
understanding  of  the  variables  and  forces  which  lead  to  the  initiation  and 
consequent  advection  of  the  dense  fog  over  KSLC. 

Direct  quantitative  comparison  between  these  studies  is  more  difficult  for 
several  reasons.  The  first  is  due  to  the  difference  in  the  J-box  used  in  each 
study.  Chilcoat  used  an  8-km  x  8-km  grid  box  summed  over  the  first  two  rj- 
levels.  In  order  to  meet  Gaussianity  requirements,  this  study  implemented  a  12- 
km  x  12-km  J-box  summed  over  the  first  two  /7-levels.  The  larger  J-box  used  in 
this  study  means  the  response  variable  includes  contributions  from  32  horizontal 
and  vertical  model  grid  points,  rather  than  the  18  used  by  Chilcoat.  This  increase 
in  J-box  size  naturally  means  the  quantified  sensitivity  response  in  this  study  will 
have  greater  magnitude. 

The  second  reason  comparison  is  difficult  is  the  change  in  DA  procedures 
that  was  undertaken  between  studies.  Assimilation  of  observations  in  this  study 
was  completed  through  ingest  of  observed  “truth”  rawindsonde  soundings  at  the 
same  time  interval  and  at  the  same  grid  spacing  as  naturally  occurred.  Chilcoat’s 
DA  process  assimilated  “truth”  rawindsonde  observations  in  columns  meant  to 
mimic  a  sparse  data  network,  contributing  100  rawinsonde  ascents  every  three 

hours.  Chilcoat’s  observations  were  randomly  located  throughout  the  domain, 
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but  fixed  in  time.  This  variation  inevitably  caused  changes  in  model  analyses 
that  would  have  resulted  in  altered  forecasts  and  thus  altered  ESA  results.  The 
specific  difference  between  the  two  systems  was  not  examined  in  this  study, 
though  its  existence  further  complicates  comparison  between  the  two  studies. 

In  addition  to  the  reasons  noted  above,  Chilcoat’s  ESA  was  performed 
over  a  different  time  frame.  Due  to  non-linearities  in  results  encountered  over  his 
original  12-hour  ESA  using  1200  UTC  23  January  analysis  and  0000  UTC  24 
January  forecast,  Chilcoat  shifted  his  12-hour  ESA  6-hours  later  in  time.  The 
ESA  results  presented  by  Chilcoat  involved  1800  UTC  23  January  analysis  and 
0600  UTC  24  January  forecast.  This  time  frame  shift  allowed  for  the  avoidance 
of  the  non-linearities  introduced  in  the  u-  and  v-  wind  components  by  the  mid¬ 
level  disturbance.  This  in  turn  allowed  for  more  linear  sensitivity  responses  than 
those  presented  in  this  study. 

B.  POTENTIAL  NETWORK  DESIGN 

The  ESA  conducted  for  this  thesis  identified  several  regions  of  sensitivity 
importance.  These  locations,  identified  in  Figure  42,  represent  locations  of 
maximum  sensitivity  or  locations  that  exhibit  dynamic  sensitivity  over  the  12-hour 
time  frame  leading  up  to  the  dense  fog  event.  It  is  believed  that  additional 
observations  placed  in  these  locations  would  provide  improved  accuracy  in  the 
prediction  of  dense  fog  events  at  KSLC.  Region  1  in  Figure  42  represents  the 
region  of  strongest  sensitivity  and  sensors  placed  here  would  be  critical  to 
observing  any  warm  air  advection  being  channeled  northward  through  the 
Central  Utah  Desert  Range.  Region  2  in  Figure  42  annotates  a  dynamic  area 
where  southerly  warm  air  advection  meets  the  influence  of  the  GSL  cold  dome. 
Additional  observations  placed  in  this  region,  the  southeast  Salt  Lake  Valley,  or 
along  the  Traverse  Mountains  to  the  south  would  help  to  define  the  warm  air 
advection  as  well  as  the  PBL  profile  in  the  southeast  GSL  Basin,  critical  for 
forecasting  dense  fog  at  KSLC.  Region  3  in  Figure  42  represents  a  low  lying 
area  northeast  of  the  GSL  which  exhibits  strong  sensitivity  when  significant 
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warm-air  advection  is  present  across  the  GSL  basin.  Though  a  single 
observation  exists  on  the  northeast  boundary  of  this  region  at  Brigham  City 
Airport  (ICAO:  KBMC),  additional  observations  across  this  plain  would  allow 
more  thorough  observation  of  warm  air  channeling  between  the  Wasatch  Range 
and  a  cold  dome  over  the  GSL.  Combined,  observations  at  these  areas  would 
help  to  characterize  the  PBL  conditions  in  the  GSL  known  to  lead  to  dense  fog 
events  at  KSLC. 
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Figure  42.  Locations  of  strong  and/or  reliable  temperature  sensitivity  regions 
that  are  considered  important  in  forecasting  qv  at  Salt  Lake  City 
International  Airport  (After  Google  Maps). 


Chilcoat  (2012)  used  u-  and  v-  wind  components  to  identify  two  regions  of 


sensitivity  he  believed  were  most  influential  to  forecasting  fog  at  KSLC  (Figure 
43).  As  discussed,  Chilcoat’s  results  represent  conclusions  from  an  ESA 
focused  6  hours  later  than  this  study.  Because  of  this  disparity  his  results  likely 


93 


represent  locations  where  dense  fog  advection  forcing  mechanisms  are 
observed,  and  not  locations  where  PBL  conditions  leading  to  dense  fog  creation 
are  focused. 


Figure  43.  Chilcoat’s  map  of  locations  of  the  strong  sensitivities  that  are 


considered  important  in  forecasting  the  amount  of  moisture  received 
at  the  SLC  International  Airport.  (From  Chilcoat  2012) 

Validation  and  confirmation  of  these  potential  network  observation 
placement  regions  would  be  a  primary  focus  of  future  work  involving  other  dense 
fog  events.  We  believe  temperature  as  ESA  initial  condition  provides  a  more 
resilient  response  to  non-linear  influence  over  longer  time  scales,  though  u-  and 
v-component  wind-based  ESA  should  also  provide  a  response  helpful  for 
complete  understanding  of  dense  fog  event  evolution. 


94 


C.  APPLICABILITY  TO  U.S.  MILITARY 

As  technology  continues  to  move  forward  providing  greater  diversity  and 
deployability  in  weather  sensor  technology,  it  is  important  that  attention  be  paid 
to  sensor  placement  which  maximizes  impact  to  forecaster  and  operator.  While 
garrison  fixed-base  observing  equipment  will  always  be  a  requirement,  the 
potential  benefit  of  the  introduction  of  additional  sensors  outside  the  base  proper 
holds  great  promise.  ESA  provides  a  means  by  which  the  impact  of  additional 
observations  can  be  determined  a  priori,  at  merely  the  cost  in  time  and  resources 
to  conduct  the  study.  As  budgetary  constraints  will  always  limit  the  amount  and 
type  of  observing  capability  available  to  weather  personnel,  this  system  provides 
a  means  of  intelligently  placing  the  limited  resources  that  are  made  available. 
The  difficulty  may  arise  when  an  ESA-determined  sensor  placement  location  is 
identified  that  provides  critical  input  to  an  NWP  system,  while  adding  limited 
additional  information  to  the  forecaster  on  the  ground.  The  decision  on  where 
and  when  to  place  a  sensor  must  be  weighed  between  the  potential 
enhancement  of  situational  awareness  provided  by  placing  a  sensor(s)  in  a 
known  data  void  region,  versus  the  potential  gain  in  forecast  accuracy  achieved 
through  placement  of  a  sensor(s)  that  primarily  benefits  the  NWP  system.  While 
NWP  will  never  provide  a  sole  source  of  informing  forecasters  in  providing 
decision-quality  information,  improving  the  data  that  it  does  provide  should  be  a 
primary  concern  to  the  Air  Force  weather  community. 

The  computing  resources  necessary  to  perform  this  ensemble  analysis 
are  significant  and  prohibit  the  implementation  of  ESA  for  data-thinning  during 
data  assimilation  or  real-time  updates  to  forecast  metrics.  The  ensemble 
computations  completed  at  ERDC  took  several  days  to  complete  and  multiple 
delays  and  issues  were  encountered.  As  computing  efficiency  increases  and 
ensemble  systems  become  more  capable  over  time,  operational  use  of  ESA  for 
these  functions  may  prove  possible.  Implementation  of  ESA  at  this  time  is  limited 
to  preconceived  sensitivity  studies  as  of  the  type  examined  in  this  study. 
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D. 


FUTURE  WORK 


This  study  showed  the  efficacy  of  using  temperature  as  a  predictor 
variable  when  qv  is  the  associated  predictand  forecast  variable.  These  results 
are  unique  to  this  mesoscale  environment  and  model  grid  scale.  Several  ideas 
are  presented  below  which  identify  means  of  further  validating  the  conclusions 
presented  in  this  paper  with  the  goal  of  expanding  the  implications  of  this 
research  to  future  and  operational  environments. 

1 .  Expand  the  research  to  include  results  for  multiple  dense  fog  events  at 
KSLC  and  across  the  GSL.  More  definitive  conclusions  regarding  the  use  of 
ESA  to  design  future  observational  networks  can  be  achieved  through 
examination  of  more  than  one  dense  fog  event.  Specific  consideration  should  be 
given  to  events  that  meet  known  climatological  norms  for  dense  fog  events.  It 
was  discovered  during  this  study  that  temperature  at  KSLC  was  >  0°C  during  the 
dense  fog  event  of  23-24  January.  This  anomaly  occurs  in  <  6%  of  all  known 
dense  fog  events  at  KSLC  over  the  past  30  years.  The  selection  of  dense  fog 
events  that  fall  within  climatological  norms  will  help  to  ensure  the  results  are 
representative  of  typical  mesoscale  forcing  in  and  around  GSL.  Study  of  other 
dense  fog  events  should  focus  on  initial  condition  PBL  temperature  and  u  and  v- 
wind  parameters,  while  also  examining  the  potential  sensitivity  response  from 
initial  condition  qvand  mu. 

2.  Evaluate  the  impact  of  multiple  additional  observations  at  the  locations 

identified  as  the  most  sensitive  grid  points.  This  study  examined  the  model 

analysis  and  forecast  change  attributable  to  the  addition  of  a  single  additional 

observation.  By  theoretically  adding  multiple  sensors  and  determining  individual 

and  collective  grouping  response,  we  can  determine  a  more  complete 

understanding  of  what  an  observational  network  design  should  resemble.  This 

evaluation  can  be  done  through  assimilating  observations  at  proposed  sites 

following  ESA  testing  procedures  for  “truth”  observations  outlined  in  section 

IV.C.2.  of  this  paper.  Following  DA,  the  analysis  difference  created  (6x)  is 

multiplied  times  the  calculated  sensitivity  ( dJ/dx )  at  that  point  to  provide  an 
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expected  forecast  change  in  qv.  It  is  not  necessary  to  rerun  the  ensemble  to 
create  a  new  expected  forecast.  This  process  is  then  repeated  at  multiple  sites 
to  test  the  design  of  a  potential  observational  network. 

3.  Place  observing  sensors  in  the  identified  sensitivity  locations  for  a  full 
winter  season.  Perform  data  denial  experiments  to  confirm  theoretical  findings 
identified  in  this  study.  This  study  utilized  synthetic  perfect  (no  error) 
observations  to  test  the  results  of  the  ESA.  Assimilation  of  real-world 
observations  introduces  additional  observation  error  which  may  affect  sensitivity 
testing  outcomes.  A  field  test  of  these  sensitivity  findings  could  help  confirm  or 
deny  the  results  found  in  this  study. 

4.  Examine  the  impact  of  alternate  observing  platforms  at  the  most 
sensitive  grid  point.  This  study  examined  how  temperature  measured  at  near¬ 
surface  level  at  points  around  the  GSL  contributed  to  formation  of  fog  at  KSLC  at 
a  later  time.  The  conclusion  from  this  study  is  that  several  fixed  sensors  or  a 
single  mobile  sensor  would  need  to  be  implemented  to  realize  the  forecast 
accuracy  improvement  suggested  by  this  ESA.  However,  it  is  unknown  exactly 
what  platform  could  provide  the  data  necessary  to  achieve  these  results.  Does  a 
fixed  sensor  that  provides  only  surface  observations  provide  enough  data  to 
sufficiently  inform  the  model  to  improve  the  forecast?  If  the  fixed  observation  is 
tower-mounted,  what  level  must  it  reach  to  gather  necessary  data?  Alternately, 
could  the  same  or  better  results  be  obtained  through  the  use  of  a  UAS  or 
aerostat-like  system?  These  systems  could  be  deployed  temporarily  based  upon 
known  synoptic  or  mesoscale  precursors  that  suggest  development  of  certain 
forecast  conditions  is  likely. 

5.  Investigate  other  potential  dense  fog  parameters  for  J  to  ensure  PBL 
moisture  budget  is  being  properly  accounted  for  and  to  maximize  the  accuracy  of 
the  fog  forecast.  A  notable  difference  (~1 .0x1 0"4  kg  kg"1  -  1 .5  kg  kg'1)  was  found 
between  observed  qv  at  KSLC  and  the  ensemble  mean  forecast  value  at 
between  1800  UTC  23  January  and  0600  UTC  24  January.  This  discrepancy 
represents  -25%  difference  between  forecast  and  observed  qv  and  may 
represent  the  cause  for  the  differential  between  actual  and  observed  values 
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obtained  during  ESA  testing.  Accurate  representation  of  the  PBL  energy  and 
moisture  budget  is  important  to  understanding  the  creation  and  dissolution  of 
dense  fog  in  the  GSL. 
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